#include <stdio.h>
#include <stdlib.h>
#include <dlfcn.h>
#include <stdio.h>
#include <iostream>
#include <vector>
#include "math.h"
#include <fstream>
#include <iostream>
#include <utility>
#include <algorithm>
#include <math.h>
#include <iostream>
#include <vector>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <iterator>
#include <vector>
#include <iostream>
#include <sstream>
#include <string>
#include <cstdlib>
#include "math.h"
#include <stdio.h>
#include <time.h>
#include <algorithm>
//#include <curses.h>
#include "TROOT.h"
#include "TFile.h"
#include <TMatrixDSym.h>
#include <TMatrixDSymEigen.h>
#include <TVectorD.h>
#include "TF1.h"
#include <TStyle.h>
#include <TSystem.h>
#include <TChain.h>
#include "TCanvas.h"
#include "TH1F.h"
#include "TLegend.h"
#include "AnalyzeSA.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "TObject.h"
#include "TClonesArray.h"
#include "TChain.h"
#include "TString.h"
#include "THStack.h"
#include "TLegend.h"
#include "TPaveText.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "Math/GenVector/VectorUtil.h"
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/Point3D.h"
#include "TH2.h"
#include "TH1.h"
#include "TH2F.h"
#include "TH1F.h"
#include "TPaveLabel.h"
#include "TPaveText.h"
#include "ExRootAnalysis/ExRootAnalysis/ExRootClasses.h"
#include "ExRootAnalysis/ExRootAnalysis/ExRootTreeReader.h"
#include "ExRootAnalysis/ExRootAnalysis/ExRootTreeWriter.h"
#include "ExRootAnalysis/ExRootAnalysis/ExRootTreeBranch.h"
#include "ExRootAnalysis/ExRootAnalysis/ExRootResult.h"
#include "ExRootAnalysis/ExRootAnalysis/ExRootUtilities.h"
/*
 #include "ExRootAnalysis/ExRootClasses.h"
 #include "ExRootAnalysis/ExRootTreeReader.h"
 #include "ExRootAnalysis/ExRootTreeWriter.h"
 #include "ExRootAnalysis/ExRootTreeBranch.h"
 #include "ExRootAnalysis/ExRootResult.h"
 #include "ExRootAnalysis/ExRootUtilities.h"
 */
#include "TApplication.h"
#include "observables.h"
#include "filters.h"
#include "jetclustering.h"
#include "hWW.h"
#include "dijet.h"
#include "utility.h"
//#include "services.h"
#include "TRandom1.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include "TStopwatch.h"
#include "TF1.h"
#include "TUnuran.h"
#include "TUnuranContDist.h"

#define normal "\033[39m"
#define	cyan "\033[36m"
#define black "\033[30m"
#define red "\033[31m"
#define green "\033[32m"
#define yellow "\033[33m"
#define blue "\033[34m"
#define magenta "\033[35m"
#define white "\033[37m"
//
//char blue[] = { 0x1b, '[', '1', ';', '3', '4', 'm', 0 };
//	char ublue[] = { 0x1b, '[', '4', ';', '3', '4', 'm', 0 };
//	char rblue[] = { 0x1b, '[', '7', ';', '3', '4', 'm', 0 };
//	char magenta[] = { 0x1b, '[', '1', ';', '3', '5', 'm', 0 };
//	char red[] = { 0x1b, '[', '1', ';', '3', '1', 'm', 0 };
//	char normal[] = { 0x1b, '[', '0', ';', '3', '9', 'm', 0 };
//	char black[] = { 0x1b, '[', '0', ';', '3', '0', 'm', 0 };
//	char green[] = { 0x1b, '[', '0', ';', '3', '2', 'm', 0 };
//	char yellow[] = { 0x1b, '[', '0', ';', '3', '3', 'm', 0 };
//	char purple[] = { 0x1b, '[', '0', ';', '3', '5', 'm', 0 };
//	char cyan[] = { 0x1b, '[', '0', ';', '3', '6', 'm', 0 };
//	char Lgray[] = { 0x1b, '[', '0', ';', '3', '7', 'm', 0 };
//	char Dgray[] = { 0x1b, '[', '0', ';', '3', '8', 'm', 0 };
//	char Bred[] = { 0x1b, '[', '1', ';', '3', '1', 'm', 0 };

using namespace std;

//extern "C" void hello(void) {
//	printf("Hello, world!\n");
//}

struct plotsrequired {
	bool sum;
	bool bg;
	bool components;
	bool normalized;
	bool diff;
};

bool check_range_request(string req_range) {
	vector<float> upperbound, lowerbound;
	bool valid_range = false;
	string override_flag;
	string override_range;
	vector<float> manual_range;
	override_flag = req_range.substr(0, 8);
	Int_t NbinFix;
	vector<float> vLBbins;

	if (strcmp(req_range.c_str(), "full") == 0) {
		valid_range = true;

	} else if (strcmp(req_range.c_str(), "optimal") == 0) {
		valid_range = true;

	} else if (strcmp(override_flag.c_str(), "override") == 0) {

		size_t posrange;
		posrange = req_range.find_first_of(",") + 1; // position of the last "," in the CLI
		override_range = req_range.substr(posrange);
		//cout<<override_range<<endl;
		//cout << "plotting " << names[aproc] << " in the manual range: ";

		float d, minimo, massimo, c, bsize;
		manual_range = get_CSVf(override_range);
		//printvector(manual_range);

		minimo = manual_range[0];
		massimo = manual_range[1];
		bsize = manual_range[2];
		c = (massimo + minimo) / 2;
		//adjust minimo and massimo according to the binsize provided by manual_range[2]
		int nofbins;
		nofbins = ceil((massimo - minimo) / bsize);
		minimo = c - nofbins * bsize / 2;
		massimo = c + nofbins * bsize / 2;
		cout << minimo << "~" << massimo << endl;
		lowerbound.push_back(minimo);
		upperbound.push_back(massimo);

		if (nofbins > 0) {
			valid_range = true;
		}

	} else if (strcmp(override_flag.c_str(), "fixedbin") == 0) {
		vector<float> vLBbins;
		size_t posrange;
		posrange = req_range.find_first_of(",") + 1; // position of the last "," in the CLI
		override_range = req_range.substr(posrange);
		//			hello();
		vLBbins = get_CSVf(override_range);
		//			printvector(vLBbins);
		NbinFix = vLBbins.size() - 1;
		lowerbound.push_back(vLBbins[0]);
		upperbound.push_back(vLBbins[vLBbins.size()]);

		if (!(vLBbins.empty())) {
			valid_range = true;
		}

	} else {
		puts("fixed wrong speicification of the range ");
		exit(1);

	}
	return valid_range;
}

//****************************************************************************
void analyse_proclist(vector<vector<string> > & processes,
		vector<double> & sigmas, vector<double> & targetcrosssections,
		vector<string> & names, vector<string> & cut_obs,
		vector<vector<int> > & cut_PID, vector<int> & cut_k,
		vector<float> & cut_lower, vector<float> & cut_upper,
		vector<int> & cut_howmany, vector<string> & cut_howmany_relation,
		vector<string> & cut_policy, vector<string> & cut_distinctpolicy,
		vector<vector<string> > & filtri_single_particles_cut,
		map<string, float> & mymap,
		map<string, float> & jetdefinition_parameters_map,
		map<string, string> & map_plot,
		vector<string> & filtri_singola_particella, vector<string> & filtri_k,
		vector<vector<string> > & filtri_k_cut) {
	//****************************************************************************

	//	char blue[] = { 0x1b, '[', '1', ';', '3', '4', 'm', 0 };
	//	char ublue[] = { 0x1b, '[', '4', ';', '3', '4', 'm', 0 };
	//	char magenta[] = { 0x1b, '[', '1', ';', '3', '5', 'm', 0 };
	//	char red[] = { 0x1b, '[', '1', ';', '3', '1', 'm', 0 };
	//	char normal[] = { 0x1b, '[', '0', ';', '3', '9', 'm', 0 };

	ifstream proclist("analysis.human");
	string commentchar;
	commentchar = "#";
	bool plot_session_open = false;
	bool proc_session_open = false;
	bool arranger_session_open = false;
	bool cut_session_open = false;
	bool analysis_session_open = false;
	bool process_open = false;
	bool jet_session_open = false;
	bool single_filter_session_open = false;
	bool k_filter_session_open = false;
	size_t found;

	vector < string > lines;
	vector < string > fields;
	string temp;
	string tempfield;
	istringstream infield;
	istringstream instream;
	vector < string > thisprocess;
	int cutcounter = 0;
	bool first = true;
	while (getline(proclist, temp)) {
		size_t index_first_not_space;
		index_first_not_space = temp.find_first_not_of(" ");

		if (temp.compare(index_first_not_space, 1, "#") != 0) {
			// since this line is not a comment it has to be divided in fields and analysed

			lines.push_back(temp);
			instream.clear();
			instream.str(temp);
			fields.clear();

			while (getline(instream, tempfield, ' ')) {
				//cout<<"tempfield is " <<tempfield<<endl;
				if (tempfield != "" && tempfield != " ") {
					fields.push_back(tempfield);
					//cout<<"filling fields with "<<tempfield<<endl;
				}
			}
			//all fields are now filled
			//the first field has to be cleaned from unwanted spaced at the end.

			fields[0] = fields[0].substr(0, fields[0].find_first_of(" "));

			//check if this is a  END ANALYSIS SECTION line
			if (fields[0].compare("*END_PLOT*") == 0) {
				plot_session_open = false;
				cout << "end of plot section reached" << endl;
			}

			if (plot_session_open == true) {

				//map<string, string>::iterator itS;
				map_plot.insert(pair<string, string>(fields[0], fields[1]));
			}

			if (fields[0].compare("*NEW_PLOT*") == 0) {
				plot_session_open = true;
				cout << "found new plot section" << endl;
			}

			//check if this is a  END ANALYSIS SECTION line
			if (fields[0].compare("*END_ANALYSIS*") == 0) {
				analysis_session_open = false;
				cout << "end of analysis section reached" << endl;
			}

			if (analysis_session_open == true) {

				map<string, float>::iterator it;
				mymap.insert(
						pair<string, float>(fields[0],
								atof(fields[1].c_str())));
			}

			if (fields[0].compare("*NEW_ANALYSIS*") == 0) {
				analysis_session_open = true;
				cout << "found new analysis section" << endl;
			}

			//check if this is a  END ARRANGER SECTION line
			if (fields[0].compare("*END_K_ARRANGE*") == 0) {
				arranger_session_open = false;
				cout << "end of arranger section reached" << endl;
			}

			if (arranger_session_open == true) {
				//arrangers.push_back(fields[0]);
			}

			if (fields[0].compare("*BEGIN_K_ARRANGE*") == 0) {
				arranger_session_open = true;
				cout << "found new arranger section" << endl;
			}

			//check if this is a  END SINGLE FILTER SECTION line
			if (fields[0].compare("*END_SINGLE_FILTERS*") == 0) {
				single_filter_session_open = false;
				cout << "end of single particle filter section reached" << endl;
			}

			if (single_filter_session_open == true) {
				filtri_singola_particella.push_back(fields[0]);
			}

			if (fields[0].compare("*BEGIN_SINGLE_FILTERS*") == 0) {
				single_filter_session_open = true;
				cout << "found new single particle filter section" << endl;
			}

			//check if this is a  END K-FILTER SECTION line
			if (fields[0].compare("*END_K_FILTERS*") == 0) {
				k_filter_session_open = false;
				cout << "end of k-filters section reached" << endl;
			}

			if (k_filter_session_open == true) {
				filtri_k.push_back(fields[0]);
			}

			if (fields[0].compare("*BEGIN_K_FILTERS*") == 0) {
				k_filter_session_open = true;
				cout << "found new kplets filter section" << endl;
			}

			if (fields[0].compare("*END_CUT_SESSION*") == 0) {
				cut_session_open = false;
				cout << "end of cut section reached" << endl;
			}
			//fill CUT variables
			if (cut_session_open == true) {
				int lastfield;
				lastfield = fields.size() - 1;
				vector<int> PIDS;
				//printvector(fields);
				cut_obs.push_back(fields[0]);

				/* */
				string inPIDstream;
				inPIDstream = fields[1];
				vector < string > PIDandFilters;
				PIDandFilters = get_USstrings(inPIDstream);
				/* */
				cut_PID.push_back(get_CSV(PIDandFilters[0]));
				filtri_single_particles_cut.push_back(
						stripelement(0, PIDandFilters));

				cut_k.push_back(atoi(fields[2].c_str()));
				//fields[3] will be the list of kplets filter
				/* */
				string inkfilters;
				inkfilters = fields[3];
				vector < string > kFilters;
				kFilters = get_USstrings(inkfilters);
				/* */
				filtri_k_cut.push_back(kFilters);

				cut_lower.push_back(atof(fields[4].c_str()));
				cut_upper.push_back(atof(fields[5].c_str()));
				// here take the request on the number of objects that must pass the cut.
				// N>k or N<k or N=k or N=-1 which means all
				cut_howmany_relation.push_back(fields[6]);
				cut_howmany.push_back(atoi(fields[7].c_str()));
				//cout<<fields[7]<<endl;
				cut_policy.push_back(fields[8]);
				cut_distinctpolicy.push_back(fields[9]);

				if ((strcmp(fields[8].c_str(), "loose") != 0)
						&& (strcmp(fields[8].c_str(), "strict") != 0)) {
					cout << fields[8]
							<< " is not a valid counting policy for a cut"
							<< endl;
					exit(1);
				}

				if ((strcmp(fields[9].c_str(), "any") != 0)
						&& (strcmp(fields[9].c_str(), "distinct") != 0)) {
					cout << fields[9]
							<< " is not a valid identity policy for a cut"
							<< endl;
					exit(1);
				}

				//printvector(cut_howmany_relation);
				string cutoutput;

				cutoutput += "new ";
				cutoutput += cut_policy[cutcounter];
				cutoutput += " cut ";
				cutoutput += float2string(cut_lower[cutcounter]);
				cutoutput += "<=";
				cutoutput += cut_obs[cutcounter];
				cutoutput += "<=";
				cutoutput += float2string(cut_upper[cutcounter]);
				cutoutput += " for ";
				cutoutput += cut_distinctpolicy[cutcounter];
				cutoutput += " N";
				cutoutput += cut_howmany_relation[cutcounter];
				cutoutput += int2string(cut_howmany[cutcounter]);
				cutoutput += " ";
				cutoutput += int2string(cut_k[cutcounter]);
				cutoutput += "-plet(s) built out of PIDs: ";
				cutoutput += vector2string(cut_PID[cutcounter]);

				cout << blue << cutoutput << normal << endl;

				//cout << ublue << "new " << cut_policy[cutcounter] << " cut " << cut_lower[cutcounter] << "<=" << cut_obs[cutcounter] << "<=" << cut_upper[cutcounter] << " for N" << cut_howmany_relation[cutcounter] << cut_howmany[cutcounter] << " " << cut_k[cutcounter] << "-plet(s) " << " built out of PIDs: ";
				//cout << magenta;
				//printvector(cut_PID[cutcounter]);
				//cout << normal;
				cutcounter++;

			}

			if (fields[0].compare("*CUT_SESSION*") == 0) {
				cut_session_open = true;
				cout << "found new cut section" << endl;

			}

			//check if this is a PROCESS SECTION line
			if (fields[0].compare("*END_PROC_SESSION*") == 0) {
				proc_session_open = false;
				cout << "end of process section reached" << endl;
			}

			if (fields[0].compare("*PROC_SESSION*") == 0) {
				proc_session_open = true;
				cout << "found new process section" << endl;
			}

			if (fields[0].compare("*ENDPROC*") == 0) {
				//cout<<"end of this process"<<endl;
				processes.push_back(thisprocess);
				thisprocess.clear();
				process_open = false;
			}

			if ((proc_session_open == true) && (process_open == true)) {
				//cout<<process_open<<process_open<<endl;
				thisprocess.push_back(fields[0]);
				//cout << fields[0] << " added to the chain" << endl;
			}

			if (proc_session_open == true
					&& fields[0].compare("*NEWPROC*") == 0) {
				process_open = true;
				double targetcrosssection = 0;
				string weightmessage;

				if (isinstring(fields[1], "pb")) {

					string inputstr = fields[1];
					inputstr.erase(
							std::remove(inputstr.begin(), inputstr.end(), 'p'),
							inputstr.end());
					inputstr.erase(
							std::remove(inputstr.begin(), inputstr.end(), 'b'),
							inputstr.end());

					targetcrosssection = atof(inputstr.c_str());
					weightmessage = " with target cross-section = " + inputstr
							+ " pb";

				} else {
					weightmessage = " with rescaling factor = " + fields[1];
				}

				sigmas.push_back(atof(fields[1].c_str()));
				targetcrosssections.push_back(targetcrosssection);
				names.push_back(fields[2]);
				cout << "new process " << red << fields[2] << normal
						<< weightmessage << " encoutered" << endl;
			}

			if (fields[0].compare("*END_JET_DEFINITION*") == 0) {
				jet_session_open = false;
				cout << "end of jet section reached" << endl;
			}

			if (jet_session_open == true) {
				float alg;
				//cout<<fields[0]<< " "<< atof(fields[1].c_str())<<endl;

				if (fields[0].compare("algorithm") == 0) {

					if (fields[1].compare("jade") == 0) {
						//						cout<<"faremo JADE"<<endl;
						alg = 0;
					}
					if (fields[1].compare("kt") == 0) {
						//						cout<<"faremo kT"<<endl;
						alg = 1;
					}

					//						cout<< "found alg"<<alg<<endl;

					jetdefinition_parameters_map.insert(
							pair<string, float>(fields[0], alg));

				} else {
					jetdefinition_parameters_map.insert(
							pair<string, float>(fields[0],
									atof(fields[1].c_str())));
				}
			}

			if (fields[0].compare("*NEW_JET_DEFINITION*") == 0) {
				jet_session_open = true;
				cout << "found new jet section" << endl;
			}

		}
	}
}

//****************************************************************************
void histo2txt(TH1 *histoE, char *txt_file) {
	//****************************************************************************
	int ns = 0;
	ns = histoE->GetNbinsX();
	Double_t xs[ns], ys[ns], exs[ns], eys[ns];
	//char* output_file="importing.txt";
	ofstream spectxt;
	spectxt.open(txt_file);

	for (Int_t i = 1; i <= ns; i++) {
		//xs[i - 1] = (histoE->GetBinCenter(i));// / log10(exp(1));
		xs[i - 1] = (histoE->GetBinLowEdge(i)); // / log10(exp(1));
		ys[i - 1] = histoE->GetBinContent(i);
		//exs[i-1]=sqrt(histoE->GetBinContent(i));
		//eys[i-1]=0.5*histoE->GetBinWidth(1);
		spectxt << xs[i - 1] << " " << ys[i - 1] << endl;
	}
	spectxt.close();
	cout << green << "Histogram output to text file: " << txt_file << normal
			<< endl;
	cout << cyan << "Import[\"" << txt_file << "\", \"Table\"]" << normal
			<< endl;

}

//****************************************************************************
void histo2txtint(TH1 *histoE, char *txt_file) {
	//****************************************************************************
	int ns = 0;
	ns = histoE->GetNbinsX();
	Double_t xs[ns], ys[ns], exs[ns], eys[ns];
	//char* output_file="importing.txt";
	ofstream spectxt;
	spectxt.open(txt_file);

	for (Int_t i = 1; i <= ns; i++) {
		xs[i - 1] = (histoE->GetBinLowEdge(i)); // / log10(exp(1));
		ys[i - 1] = histoE->GetBinContent(i);
		//exs[i-1]=sqrt(histoE->GetBinContent(i));
		//eys[i-1]=0.5*histoE->GetBinWidth(1);
		if (ys[i - 1] != 0) {
			double localX = xs[i - 1];
			spectxt << any2closestinteger(localX) << " " << ys[i - 1] << endl;
		}
	}
	spectxt.close();

	cout << green << "Histogram output to text file: " << txt_file << normal
			<< endl;

}

//****************************************************************************
void histo2nbint(TH1 *histoE, char *txt_file) {
	//****************************************************************************
	int ns = 0;
	bool notfirstitem = false;
	ns = histoE->GetNbinsX();
	Double_t xs[ns], ys[ns], exs[ns], eys[ns];
	//char* output_file="importing.txt";
	ofstream spectxt;
	spectxt.open(txt_file);
	spectxt << "{";
	for (Int_t i = 1; i <= ns; i++) {
		xs[i - 1] = (histoE->GetBinLowEdge(i)); // / log10(exp(1));
		ys[i - 1] = histoE->GetBinContent(i);
		//exs[i-1]=sqrt(histoE->GetBinContent(i));
		//eys[i-1]=0.5*histoE->GetBinWidth(1);
		if (ys[i - 1] != 0) {
			double localX = xs[i - 1];
			if (notfirstitem) {
				spectxt << ",";
			}
			spectxt << "{" << any2closestinteger(localX) << "," << ys[i - 1]
					<< "}";
			notfirstitem = true;
		}
	}
	spectxt << "}" << endl;
	spectxt.close();
}

//****************************************************************************
void histo2nb(TH1 *histoE, char *txt_file) {
	//****************************************************************************
	int ns = 0;
	bool notfirstitem = false;
	ns = histoE->GetNbinsX();
	Double_t xs[ns], ys[ns], exs[ns], eys[ns];
	//char* output_file="importing.txt";
	ofstream spectxt;
	spectxt.open(txt_file);
	spectxt << "{";
	for (Int_t i = 1; i <= ns; i++) {
		xs[i - 1] = (histoE->GetBinLowEdge(i)); // / log10(exp(1));
		ys[i - 1] = histoE->GetBinContent(i);
		//exs[i-1]=sqrt(histoE->GetBinContent(i));
		//eys[i-1]=0.5*histoE->GetBinWidth(1);

		if (notfirstitem) {
			spectxt << ",";
		}
		spectxt << "{" << xs[i - 1] << "," << ys[i - 1] << "}";
		notfirstitem = true;

	}
	spectxt << "}" << endl;
	spectxt.close();
}

//*************************************************
void findclassmembers(vector<int> classpid, Int_t nmax,
		vector<TRootLHEFParticle*> vecpart, vector<int>& classmembers,
		int minstatus, int maxstatus) {
	//*************************************************
	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {
		int PID;
		PID = vecpart[jcount]->PID; //each time changes
		if (isinclass(classpid, PID)) {
			int status;
			status = vecpart[jcount]->Status;
			if (isinrange(status, minstatus, maxstatus) == 1) {
				classmembers.push_back(jcount);
			}
		}
		/*if ( !isinclass(classpid,PID)) {
		 cout<<"no interesting particle found"<<endl;;
		 }*/

	}
}

//****************************************************************************
// TAKING THE classpid VECTOR AS INPUT:
// FIRST)  FIND ALL THE PARTICLES WITH THAT PID (findclassmembers)
// SECOND) FORM AL THE POSSIBLE SUBSETS OF PARTICLES (findkplets)
//
//****************************************************************************
void fill_allkplets_for_these_PIDS(Int_t nmax, vector<int> classpid,
		vector<TRootLHEFParticle*> part,
		vector<vector<vector<int> > > & allkplets, vector<int> & classmembers,
		int minstatus, int maxstatus) {
	//****************************************************************************

	//fill the vector classmembers searching the PIDs
	findclassmembers(classpid, nmax, part, classmembers, minstatus, maxstatus);
	//cout<<"ho trovato "<<nmembers<<" membri per questa classe"<<endl;
	int nmembers;
	nmembers = classmembers.size();
	if (nmembers > 0) {
		//cout<<"appena inizializzato "<<allkplets.empty()<<endl;

		//trovo tutti i multipletti e li scrivo in allkplets
		for (int kplet = 1; kplet <= nmembers; kplet++) {

			puts("sto per fare findkplets");

			//cout<<"trovo le "<<kplet<<"tuple in "<<nmembers<<" particelle"<<endl;
			findkplets(allkplets, kplet, nmembers);
			//cout<<"trovate"<<endl;

			puts("ho fatto findkplets");

		}
	}

}

//****************************************************************************
void fill_allkplets_for_this_list(vector<vector<vector<int> > > & allkplets,
		vector<vector<int> > & classmembers, int kpletin) {
	//****************************************************************************

	int nmembers;
	nmembers = classmembers.size();

	int k;
	k = kpletin;

	if (nmembers > 0) {
		//cout<<"appena inizializzato "<<allkplets.empty()<<endl;

		/*********************************************************************************/
		// se volessi limitare l'esplosione della combinatorica dei casi in cui Newton[k,n]~100 lo dovrei mettere qui
		/*********************************************************************************/

		//trovo tutti i multipletti e li scrivo in allkplets
		//for (int kplet = 1; kplet <= nmembers; kplet++) {
		/* eh no bello mio! faccio un solo k se no e` inutile */

		//cout << "sto per fare findkplets per k=" << kplet << endl;
		//cout<<"trovo le "<<kplet<<"tuple in "<<nmembers<<" particelle"<<endl;
		if (kpletin == 0) {
			k = classmembers.size();
		}

		//		clock_t start, end;
		//		double dif;
		//		start = clock();
		findkplets(allkplets, k, nmembers);
		//		end = clock();
		//		dif = end-start;
		//		cout << "tempo: " << end<< " "<< start  << " per #=" << k << endl;
		//cout<<"trovate"<<endl;
		//puts("ho fatto findkplets");

		//}
	}

}

////****************************************************************************
//void memberstoevaluate(vector<vector<int> > & toevaluate,
//		vector<vector<int> > classmembers,
//		vector<vector<vector<int> > > allkplets, int k, int riga) {
//	//****************************************************************************
//	int colonna = 0;
//	while (colonna < allkplets[k - 1][riga].size()) {
//		toevaluate.push_back(classmembers[allkplets[k - 1][riga][colonna] - 1]);
//		colonna++;
//	}
//}//EndOf memberstoevaluate
//
//
////****************************************************************************
//void memberstoevaluatesimple(vector<vector<int> > & toevaluate, vector<vector<
//		int> > classmembers, vector<vector<int> > allkplets, int riga) {
//	//****************************************************************************
//	int colonna = 0;
//	while (colonna < allkplets[riga].size()) {
//		toevaluate.push_back(classmembers[allkplets[riga][colonna] - 1]);
//		colonna++;
//	}
//}//EndOf memberstoevaluate

/***********************************/
/***********************************/
//typedef long double (*FUN)(vector<int> , TLorentzVector[], TRootLHEFParticle*[]);
//vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part
typedef long double (*FUN)(vector<vector<int> >, vector<TLorentzVector>,
		vector<TRootLHEFParticle*>);
typedef vector<vector<int> > (*FIL)(vector<vector<int> >, vector<vector<int> >,
		vector<TLorentzVector>, vector<TRootLHEFParticle*>, string);

// the type for the pointers to functions of such arguments
/***********************************/
/***********************************/

////*************************************************
////vector<float> computeOBSonthekplet(TRootLHEFParticle *LHEparticles[], TLorentzVector vec_part[], vector<vector<vector<int> > > allkplets, int k, vector<int> classmembers, FUN OBS) {
//vector<float> computeOBSonthekplet(vector<TRootLHEFParticle*> LHEparticles,
//		vector<TLorentzVector> vec_part,
//		vector<vector<vector<int> > > allkplets, int k,
//		vector<vector<int> > classmembers, FUN OBS) {
//	//*************************************************
//
//	vector<float> results;
//	int riga = 0;
//	//cout<<"numero di "<<k<<"-etti="<<allkplets[k-1].size()<<endl;
//	//cout<<allkplets[k-1].empty()<<endl;
//
//
//	while (riga < allkplets[0].size()) {
//
//		vector<vector<int> > toevaluate;
//		//find the current members to evaluate
//		// if it was asked for 0 members do it for all the elements in the class
//		// this is typically what you want for the missing energy
//
//		memberstoevaluate(toevaluate, classmembers, allkplets, 1, riga);
//
//		//results.push_back(OBS(toevaluate, vec_part));
//		float value;
//		value = OBS(toevaluate, vec_part, LHEparticles);
//		//value.push_back(w);
//		results.push_back(value);
//		riga++;
//
//	}
//	return results;
//}//EndOf compute()
////*************************************************
////vector<float> computeOBSonthekplet(TRootLHEFParticle *LHEparticles[], TLorentzVector vec_part[], vector<vector<vector<int> > > allkplets, int k, vector<int> classmembers, FUN OBS) {
//vector<float> computeOBSonthekpletsimple(
//		vector<TRootLHEFParticle*> LHEparticles,
//		vector<TLorentzVector> vec_part, vector<vector<int> > allkplets,
//		vector<vector<int> > classmembers, FUN OBS) {
//	//*************************************************
//	/* return a vector of observable's values from a vector of kplets
//	 * this applyies to the case of simple particles when the kplets is a 1-plet
//	 * it always needs the translation done through classmembers
//	 * */
//
//	vector<float> results;
//	int riga = 0;
//	//cout<<"numero di "<<k<<"-etti="<<allkplets[k-1].size()<<endl;
//	//cout<<allkplets[k-1].empty()<<endl;
//
//
//	while (riga < allkplets.size()) {
//
//		vector<vector<int> > toevaluate;
//		//find the current members to evaluate
//		// if it was asked for 0 members do it for all the elements in the class
//		// this is typically what you want for the missing energy
//
//		/* translates each k-plet in the cardinal notation 1,2,3, in lines in the event, Event Notation*/
//		memberstoevaluatesimple(toevaluate, classmembers, allkplets, riga);
//
//		//results.push_back(OBS(toevaluate, vec_part));
//		float value;
//		value = OBS(toevaluate, vec_part, LHEparticles);
//		//value.push_back(w);
//		results.push_back(value);
//		riga++;
//
//	}
//	return results;
//}//EndOf compute()
//*************************************************
int ID_highest_histogram(vector<TH1*> broken) {
	//*************************************************
	int highest_histo_id = 0;
	float max_in_histos = 0;

	for (int khisto = 0; khisto < broken.size(); khisto++) {
		float max_this_histo;
		max_this_histo = broken[khisto]->GetBinContent(
				broken[khisto]->GetMaximumBin());

		if (max_this_histo > max_in_histos) {
			max_in_histos = max_this_histo;
			highest_histo_id = khisto;
		}
	}
	return highest_histo_id;
}

//*************************************************
bool pass_the_request(int howmany, string cut_howmany_relation, int cut_howmany,
		vector<vector<int> > kpletsthatpassthecut, string idpolicy) {
	//****************************************************************************
	bool check_distinct;
	bool enough;

	/* distinct can be used only with > and =*/

	if (!strcmp(idpolicy.c_str(), "any")) {

		check_distinct = false;

	} else if (!strcmp(idpolicy.c_str(), "distinct")) {
		check_distinct = true;

	}

	if (!strcmp(cut_howmany_relation.c_str(), ">")) {
		//		cout << "how, cut_ " << howmany << " " << cut_howmany << endl;
		//		puts(">");
		if (howmany > cut_howmany) {
			enough = true;
		} else if (howmany <= cut_howmany) {
			enough = false;
		}

		//		cout << "enough=" << enough << endl;

		if (enough == true) {
			/* check wheter they are distinct when they have to*/
			if (check_distinct == true) {
				bool is_distinct;

				// cout << "k-plets that pass the cut" << endl;
				//printmatrix(kpletsthatpassthecut);
				is_distinct = CheckOneDistinct(kpletsthatpassthecut,
						cut_howmany);

				//		cout << "esame di distinct " << is_distinct << endl;
				enough = is_distinct;
				return is_distinct;

			}
			return enough;
		}

		//		cout << "enough=" << enough << endl;
		return enough;

	} else if (!strcmp(cut_howmany_relation.c_str(), "<")) {
		//		puts("<");
		if (howmany < cut_howmany) {
			return true;
		} else if (howmany >= cut_howmany) {
			return false;
		}

	} else if (!strcmp(cut_howmany_relation.c_str(), "=")) {
		//		puts("=");

		if (check_distinct == false) {

			if (howmany == cut_howmany) {
				return true;
			} else if (howmany != cut_howmany) {
				return false;
			}
		} else if (check_distinct == true) {

			bool is_distinct_n_plus_one, is_distinct_n;

			// cout << "k-plets that pass the cut" << endl;
			//printmatrix(kpletsthatpassthecut);

			if (howmany >= cut_howmany + 1) {

				is_distinct_n_plus_one = CheckOneDistinct(kpletsthatpassthecut,
						cut_howmany);
			} else {
				is_distinct_n_plus_one = false;
			}
			if (howmany >= cut_howmany) {
				is_distinct_n = CheckOneDistinct(kpletsthatpassthecut,
						cut_howmany - 1);
			} else {
				is_distinct_n = false;
			}

			//		cout << "esame di distinct " << is_distinct << endl;
			if (is_distinct_n == true && is_distinct_n_plus_one == false) {
				return true;
				//	cout << "= distinct" << endl;
			} else if (is_distinct_n == true
					&& is_distinct_n_plus_one == true) {
				//	cout << "= and > distinct" << endl;
				return false;
			} else {
				//	cout << "too few" << endl;
				return false;
			}
		}

	} else {
		cout << "relation " << cut_howmany_relation << " is not valid";
		return false;
	}

}

void PrintCanvasMultiFormat(const char file[], TCanvas * c1, string extension,
		TH1 *histoE) {
	//*************************************************

	/* check if the histogram is empty otherwise the plot will break the program*/

	if (histoE->GetEntries() != 0) {

		vector < string > extensions;
		extensions = getCSstrings(extension);
		//cout << extensions.size();
		for (int iext = 0; iext < extensions.size(); iext++) {
			char tempfile[1000] = "./output/";

			string extension;
			size_t pos;

			pos = extensions[iext].find_last_of("."); // position of the last "." in the CLI
			extension = extensions[iext].substr(pos);

			strcat(tempfile, file);
			strcat(tempfile, extensions[iext].c_str());

			if (extension.compare(".txt") == 0) {

				histo2txt(histoE, tempfile);

			} else if (extension.compare(".txtint") == 0) {

				histo2txtint(histoE, tempfile);

			} else if (extension.compare(".nbint") == 0) {

				histo2nbint(histoE, tempfile);

			} else if (extension.compare(".nb") == 0) {

				histo2nb(histoE, tempfile);

			} else if (extension.compare(".value") == 0) {
// do nothing this has already been done at the moment of filling the array of the values
			} else {
				histoE->SetMinimum(0.0);
				c1->Print(tempfile);
			}
		}
	} else {
		cout << red << "WARNING: the histogram for " << file
				<< "was empty and therefore no plot has been produced" << endl;

	}
}

//*********************************************************************
//       *************      MULTI PANEL CANVAS     *************
TCanvas * histo2canvas(vector<TH1*> plot_components, string CLI_title,
		string plot_title, vector<string> x_label, int logy,
		vector<string> names) {
	//*********************************************************************

	TCanvas *c22 = new TCanvas("c22", "Histogram Drawing Options", 200, 10, 700,
			900);
//	TPad *pad1 = new TPad("pad1", "The pad with the function", 0.03, 0.62, 0.40,
//			0.92, 21);
//	TPad *pad2 = new TPad("pad2", "The pad with the histogram", 0.41, 0.62,
//			0.98, 0.92, 21);
	TPad *pad1 = new TPad("pad1", "The pad with the function", 0.03, 0.75, 0.97,
			0.92, 21);
	TPad *pad2 = new TPad("pad2", "The pad with the histogram", 0.03, 0.58,
			0.97, 0.75, 21);
	TPad *pad3 = new TPad("pad3", "The pad with the histogram", 0.03, 0.02,
			0.97, 0.57, 21);
	pad1->Draw();
	pad2->Draw();
	pad3->Draw();
	// Draw a global picture title
	TPaveLabel *title = new TPaveLabel(0.1, 0.94, 0.9, 0.98, CLI_title.c_str());
	title->SetFillColor(kWhite);
	title->SetTextFont(52);
	title->Draw();

	// Draw histogram hpx in first pad with the default option.
	pad3->cd();
	//	pad1->GetFrame()->SetFillColor(18);

	pad3->SetBorderSize(2);
	pad3->SetFrameFillColor(kWhite);
	pad3->SetFillColor(kWhite);
	pad3->SetBorderMode(28);
	pad3->SetGrid();

	if (logy == 1) {
		gPad->SetLogy(1);
	}

	//	plot_componentsname += "_components";
	string plot_componentsname = plot_title;

	int highest_histo_id;
	highest_histo_id = ID_highest_histogram(plot_components);
	cout << "the highest bin belong to the process " << names[highest_histo_id]
			<< endl;

	plot_components[highest_histo_id]->SetLineColor(highest_histo_id + 2);
	plot_components[highest_histo_id]->SetTitle(plot_componentsname.c_str());
	plot_components[highest_histo_id]->Draw("");

	for (int cproc = 0; cproc < plot_components.size(); cproc++) {
		//cout<<cproc<<endl;
		plot_components[cproc]->SetLineColor(cproc + 2);

		if (cproc != highest_histo_id) {
			plot_components[cproc]->Draw("same");
		}
	}

	pad1->cd();

	pad1->SetBorderSize(2);
	pad1->SetFrameFillColor(kWhite);
	pad1->SetFillColor(kWhite);
	pad1->SetBorderMode(28);

	int max_length_name = 0;
	for (int cproc = 0; cproc < plot_components.size(); cproc++) {
		if (names[cproc].size() > max_length_name) {
			max_length_name = names[cproc].size();
		}
	}

	TLegend* leg;
	leg = new TLegend(0.01, 0.99 - 0.08 * plot_components.size(), 0.99, 0.99);
	//	leg->SetTextSize(0.04);

	for (int cproc = 0; cproc < plot_components.size(); cproc++) {
		leg->AddEntry(plot_components[cproc], names[cproc].c_str(), "l");
	}

	//leg->SetHeader("The Legend Title");

	leg->Draw();

	pad2->cd();
	TPaveLabel *label2 = new TPaveLabel(0.05, 0.9, 0.92, 0.99, "cuts");

	label2->SetFillColor(kWhite);
	label2->Draw();
	pad2->SetBorderSize(2);
	pad2->SetFrameFillColor(kWhite);
	pad2->SetFillColor(kWhite);
	pad2->SetBorderMode(28);

	TPaveText * pt = new TPaveText(0.05, 0.01, 0.92, 0.90, "NDC"); // NDC sets coords
	// relative to pad dimensions
	pt->SetFillColor(0); // text is black on white
	pt->SetTextFont(40);
	pt->SetTextSize(0.05);
	pt->SetTextAlign(12);
	for (int cut_k = 0; cut_k < x_label.size(); cut_k++) {
		pt->AddText(x_label[cut_k].c_str());
	}
	pt->Draw(); //to draw your text object

	c22->Update();

	return c22;
}

//*************************************************
//    multiintelliplothisto2eps(globalresult,0,argv, names);
//[process#][event#][kplet#]
void multiintelliplothisto2eps(
		vector<vector<vector<float> > > globalresult,
		int logy,
		char **argv,
		vector<string> names,
		vector<float> weights,
		vector<vector<float> > globalweights,
		vector<vector<float> > global_weights_of_the_selected_events_of_this_process,
		string req_range, string extension, vector<string> x_label, int distro,
		string plot_title, string filenamesuffix,
		struct plotsrequired whichplots) {
	//*************************************************
	//cout<<"inizio il plot"<<endl;
	vector<float> upperbound, lowerbound;
	vector<int> points;
	vector<int> bins;

	//	char file[100] = "";
	//	strcat(file, argv[1]);
	//	strcat(file, "_");
	//	strcat(file, names[aproc].c_str());

	cout << "the weight vector has dimension " << weights.size() << endl;
	for (int ksample = 0; ksample < globalweights.size(); ksample++) {
		cout << "the component " << ksample
				<< " of globalweight vector is made of  "
				<< globalweights[ksample].size() << " elements" << endl;
	}

	for (int aproc = 0; aproc < globalresult.size(); aproc++) { //loop on the processes

		// *** HERE I MAKE THE HISTOGRAM FOR A SINGLE PROCESS ***

		TH1 *invmasssq;
		// sintax of the override mode: override,300,500,10   where numbers are xmin,xmax,binsize

		string override_flag;
		string override_range;
		vector<float> manual_range;
		override_flag = req_range.substr(0, 8);
		Int_t NbinFix;
		vector<float> vLBbins;

		if (strcmp(req_range.c_str(), "full") == 0) {
			cout << "plotting " << names[aproc] << " in the " << req_range
					<< " range";
			//hello();

			float d, minimo, massimo, c;
			//printmatrix(globalresult[aproc]);

			minimo = minmatrix(globalresult[aproc]);

			massimo = maxmatrix(globalresult[aproc]);

			d = massimo - minimo;
			c = (massimo + minimo) / 2;
			massimo = c + (1.01) * d / 2;
			minimo = c - (1.01) * d / 2;
			lowerbound.push_back(minimo);
			upperbound.push_back(massimo);
			//cout << "lowerbounds ";
			//printvector(lowerbound);

		} else if (strcmp(req_range.c_str(), "optimal") == 0) {
			cout << "plotting " << names[aproc] << " in the " << req_range
					<< " range. ";

			if (maxmatrix(globalresult[aproc]) - minmatrix(globalresult[aproc])
					> 4 * variancematrix(globalresult[aproc])) {

				//lowerbound.pop_back();
				lowerbound.push_back(
						avgmatrix(globalresult[aproc])
								- 2 * variancematrix(globalresult[aproc]));
				//upperbound.pop_back();
				upperbound.push_back(
						avgmatrix(globalresult[aproc])
								+ 2 * variancematrix(globalresult[aproc]));
			} else {
				float d, minimo, massimo, c;

				minimo = minmatrix(globalresult[aproc]);
				massimo = maxmatrix(globalresult[aproc]);
				//cout << endl << endl;
				//cout << minimo << endl;
				//cout << massimo << endl;
				d = massimo - minimo;
				c = (massimo + minimo) / 2;
				//cout << d << endl;
				//cout << c << endl;
				massimo = c + (1.01) * d / 2;
				minimo = c - (1.01) * d / 2;
				//cout << "alla fine: " << endl << minimo << endl;
				//cout << massimo << endl;
				lowerbound.push_back(minimo);
				upperbound.push_back(massimo);
				//cout << "lowerbounds ";
				//printvector(lowerbound);
			}

		} else if (strcmp(override_flag.c_str(), "override") == 0) {
			cout << "plot range of " << names[aproc] << " is set manually to ";

			size_t posrange;
			posrange = req_range.find_first_of(",") + 1; // position of the last "," in the CLI
			override_range = req_range.substr(posrange);
			//cout<<override_range<<endl;
			//cout << "plotting " << names[aproc] << " in the manual range: ";

			float d, minimo, massimo, c, bsize;
			manual_range = get_CSVf(override_range);
			//printvector(manual_range);

			minimo = manual_range[0];
			massimo = manual_range[1];
			bsize = manual_range[2];
			c = (massimo + minimo) / 2;
			//adjust minimo and massimo according to the binsize provided by manual_range[2]
			int nofbins;
			nofbins = ceil((massimo - minimo) / bsize);
			minimo = c - nofbins * bsize / 2;
			massimo = c + nofbins * bsize / 2;
			cout << minimo << "~" << massimo << endl;
			lowerbound.push_back(minimo);
			upperbound.push_back(massimo);
		} else if (strcmp(override_flag.c_str(), "fixedbin") == 0) {
			cout << "plot range of " << names[aproc] << " is set manually to ";

			size_t posrange;
			posrange = req_range.find_first_of(",") + 1; // position of the last "," in the CLI
			override_range = req_range.substr(posrange);
			//			hello();
			vLBbins = get_CSVf(override_range);
			//			printvector(vLBbins);
			NbinFix = vLBbins.size() - 1;
			lowerbound.push_back(vLBbins[0]);
			upperbound.push_back(vLBbins[vLBbins.size()]);
		} else {
			puts("fixed wrong speicification of the range ");
			exit(1);

		}

		points.push_back(
				elementsofamatrix(globalresult[aproc], lowerbound[aproc],
						upperbound[aproc]));

		//printvector(points);
		//with more than 10K points use sqrt(points)
		//from 1K to 10K use points/75
		//from 0 to 1K use 15 bins

		if (strcmp(override_flag.c_str(), "override") == 0) {
			bins.push_back(
					ceil(
							(manual_range[1] - manual_range[0])
									/ manual_range[2]));
		} else if (strcmp(override_flag.c_str(), "fixedbin") == 0) {

			bins.push_back(NbinFix);

		} else {

			if (points[aproc] < 1000) {
				bins.push_back(15);
			}
			if (points[aproc] >= 1000 && points[aproc] <= 10000) {
				int nbin = (int) points[aproc] / 75;
				bins.push_back(nbin);
				//cout<<bins.size()<<endl;
			}
			if (points[aproc] > 10000) {
				int nbin = (int) sqrt(points[aproc]);
				bins.push_back(nbin);
			}
		}

		int Ndiv = bins[aproc];

		cout << ": " << lowerbound[aproc] << " ~ " << upperbound[aproc] << endl;
		float binw = (upperbound[aproc] - lowerbound[aproc]) / Ndiv;

		string histogram_name;
		string histogram_title; // upto4fields "Histogram title;X Axis;Y Axis;Z Axis"

		histogram_name += argv[1];
		histogram_name += "_";
		histogram_name += names[aproc];
		histogram_name += filenamesuffix;

		histogram_title += argv[1];
		histogram_title += " ";
		histogram_title += names[aproc];
		histogram_title += ";";
		//histogram_title += x_label;
		histogram_title += " ";
		histogram_title += ";N/";
		histogram_title += float2string(binw);
		histogram_title += " GeV;";

		// declare the histogram
		const char * mainhistoname = histogram_name.c_str();
		const char * mainhistotitle = histogram_title.c_str();
		if (strcmp(override_flag.c_str(), "fixedbin") == 0) {

			float binsLB[vLBbins.size()];

for			(int k = 0; k < vLBbins.size(); k++) {
				binsLB[k] = vLBbins[k];
				//				printvector(vLBbins);
			}

			//			cout<<"lower edge of the bins"<<binsLB[0]<<" "<<binsLB[1]<<" "<<binsLB[2]<<endl;

			//			float lb[3];
			//			lb[0]=0;
			//			lb[1]=1;
			//			lb[2]=2;

			invmasssq
			= new TH1F(mainhistoname, mainhistotitle, NbinFix, binsLB);
		} else {
			invmasssq = new TH1F(mainhistoname, mainhistotitle, Ndiv,
					lowerbound[aproc], upperbound[aproc]);
		}

		cout << " the set of events with the selected properties is made of "
				<< globalresult[aproc].size() << " weighted events" << endl;

		cout
				<< " the set of weights corresponding to the selected events is made of "
				<< global_weights_of_the_selected_events_of_this_process[aproc].size()
				<< " elements" << endl;

		for (int ievent = 0; ievent < globalresult[aproc].size(); ievent++) {

			for (int kplet = 0; kplet < globalresult[aproc][ievent].size();
					kplet++) {
				float value;
				value = globalresult[aproc][ievent][kplet];
				float w =
						weights[aproc]
								* global_weights_of_the_selected_events_of_this_process[aproc][ievent];
				// global_weights_of_the_selected_events_of_this_process[aproc][ievent] is the weight of this event
				// weights is lumi * 1000 / Nsamples
				invmasssq->Fill(value, w);
			}
		}

		if (globalresult[aproc].size() > 0) {
			gStyle->SetCanvasColor(0);
			gStyle->SetTitleColor(0);
			gStyle->SetCanvasColor(0);
			//gStyle->SetPaperSize(26, 20);

			//		c1->Range(27.5, -13.2563, 152.5, 119.306);
			TCanvas *c1 = new TCanvas(mainhistoname, mainhistoname, 1);

			c1->Divide(1, 1);
			c1->cd(1);
			c1->SetBorderSize(2);
			c1->SetFrameFillColor(kWhite);
			c1->SetFillColor(kWhite);
			c1->SetBorderMode(28);
			c1->SetGrid();

			if (logy == 1) {
				gPad->SetLogy(1);
			}

			if (distro == 1) {
				// rewight each bin by its width

				int nbinsh;
				nbinsh = invmasssq->GetNbinsX();

				//GetBinLowEdge(Int_t bin)
				for (int binth = 1; binth <= nbinsh; binth++) {
					float binwidth;
					binwidth = invmasssq->GetBinWidth(binth);
					cout << "delta-i " << binwidth << endl;
					float newbincontent = (invmasssq->GetBinContent(binth))
							/ binwidth;
					invmasssq->SetBinContent(binth, newbincontent);
				}
			}
			invmasssq->SetLineColor(aproc + 2);
			invmasssq->Draw("");

			//break extensions into CSV and then print for each value
			//check if is a .txt or .dat out put and call the relative function
			PrintCanvasMultiFormat(histogram_name.c_str(), c1, extension,
					invmasssq);

			// ****** IMPLEMENTATION OF THE OUTPUT OF ALL THE COMPUTED VALUES *********
			// void ValuesToFile (  vector<vectot<vector<float >  > > globalresult, int aproc, string extension, string histogram_name  )

			vector < string > extensions;
			extensions = getCSstrings(extension);

			for (int iext = 0; iext < extensions.size(); iext++) {

				string ithextension;
				size_t pos;

				pos = extensions[iext].find_last_of("."); // position of the last "." in the CLI
				ithextension = extensions[iext].substr(pos);

				if (ithextension.compare(".value") == 0) {

					ofstream valuestxt;
					string valuesname = "./output/" + histogram_name
							+ ".value.txt";
					const char * values_filename = valuesname.c_str();
					valuestxt.open(values_filename);

					for (int ievent = 0; ievent < globalresult[aproc].size();
							ievent++) {

						for (int kplet = 0;
								kplet < globalresult[aproc][ievent].size();
								kplet++) {
							float value;
							value = globalresult[aproc][ievent][kplet];

							valuestxt << value << endl;

						}

					}
					valuestxt.close();
				} //end of the check of the extension .value

			}

			//cout<<"Number of "<<plotpid[zcount-6]<<": "<<LoghistoE->GetEntries()<<endl;
		}

	} //cycle on the processes

	// ****************************************************
	//********** SUMMARY HISTOGRAMS ************************
	// ****************************************************

	float absolutelowerbound = minvector(lowerbound);
	float absoluteupperbound = maxvector(upperbound);
	cout << absolutelowerbound << " " << absoluteupperbound << endl;
	int binnumber = minvector(bins);
	float binsize_summary_histos = (absoluteupperbound - absolutelowerbound)
			/ binnumber;

	//********** DECLARE HISTO OF THE SUM OF THE PIECES ************************

	string sum_name;
	string sum_title;
	string bg_name;
	string bg_title;

	sum_name += argv[1];
	sum_name += " ";
	sum_name += "sum";
	sum_name += filenamesuffix;

	sum_title += argv[1];
	sum_title += " ";
	sum_title += "sum";
	sum_title += ";;N/";
	sum_title += float2string(binsize_summary_histos);
	sum_title += " GeV;";

	bg_name += argv[1];
	bg_name += " ";
	bg_name += "BG";
	bg_name += filenamesuffix;

	bg_title += argv[1];
	bg_title += " ";
	bg_title += "BG";
	bg_title += ";;N/";
	bg_title += float2string(binsize_summary_histos);
	bg_title += " GeV;";

	TH1* sum;
	sum = new TH1F(sum_name.c_str(), sum_title.c_str(), binnumber,
			absolutelowerbound, absoluteupperbound);
	TH1* bg;
	bg = new TH1F(bg_name.c_str(), bg_title.c_str(), binnumber,
			absolutelowerbound, absoluteupperbound);

	//********** DECLARE HISTOS OF THE BROKEN PIECES ************************

	vector<TH1*> plot_components;
	vector<TH1*> plot_components_unit;
	string plot_componentsname;
	plot_componentsname += argv[1];
	string plot_componentsunitname;
	plot_componentsunitname += argv[1];
	plot_componentsunitname += "_normalized";

	for (int dproc = 0; dproc < globalresult.size(); dproc++) {

		string plot_componentsname_ith;
		string plot_componentstitle_ith;
		plot_componentsname_ith += argv[1];
		plot_componentsname_ith += "_";
		plot_componentsname_ith += int2string(dproc);

		plot_componentstitle_ith += argv[1];
		plot_componentstitle_ith += " ";
		plot_componentstitle_ith += names[dproc];
		plot_componentstitle_ith += ";;N/";
		plot_componentstitle_ith += float2string(binsize_summary_histos);
		plot_componentstitle_ith += " GeV;";

		string plot_componentsunitname_ith;
		string plot_componentsunittitle_ith;
		plot_componentsunitname_ith = plot_componentsname_ith + "_normalized";
		plot_componentsunittitle_ith = plot_componentstitle_ith;

		plot_components.push_back(
				new TH1F(plot_componentsname_ith.c_str(),
						plot_componentstitle_ith.c_str(), binnumber,
						absolutelowerbound, absoluteupperbound));
		plot_components_unit.push_back(
				new TH1F(plot_componentsunitname_ith.c_str(),
						plot_componentsunittitle_ith.c_str(), binnumber,
						absolutelowerbound, absoluteupperbound));
	}

	// ********* FILL HISTOS OF THE SUM ************************
	// AND
	// ********* FILL THE HISTO OF THE BROKEN PIECES ************************

	for (int bproc = 0; bproc < globalresult.size(); bproc++) {
		for (int ievent = 0; ievent < globalresult[bproc].size(); ievent++) {
			for (int kplet = 0; kplet < globalresult[bproc][ievent].size();
					kplet++) {
				float value;
				value = globalresult[bproc][ievent][kplet];

				float w =
						weights[bproc]
								* global_weights_of_the_selected_events_of_this_process[bproc][ievent];
				//global_weights_of_the_selected_events_of_this_process[aproc][ievent] is the weight of this event
				// weights is lumi * 1000 / Nsamples
				sum->Fill(value, w);
				plot_components[bproc]->Fill(value, w);
				plot_components_unit[bproc]->Fill(value, w);

				size_t found;
				found = names[bproc].find("signal");
				if (found == string::npos) {
					//cout<<names[bproc] << " is a background"<<endl;
					//if the name do not contains "signal" add it to the BG
					bg->Fill(value, w);
				}

			}
		}
		//cout<<"peso di "<<names[bproc]<<" "<<weights[bproc]<<endl;
	}

	int IDhistoA, IDhistoB;
	for (int bproc = 0; bproc < plot_components.size(); bproc++) {

		size_t found;

		found = names[bproc].find("diffA");
		if (!(found == string::npos)) {
			IDhistoA = bproc;
		}

		found = names[bproc].find("diffB");
		if (!(found == string::npos)) {
			IDhistoB = bproc;

		}
	}

	if (whichplots.diff == true) {
		string diff_name, diff_title;

		diff_name += argv[1];
		diff_name += " ";
		diff_name += "deviation";
		diff_name += filenamesuffix;

		diff_title += argv[1];
		diff_title += " ";
		diff_title += "deviation";
		diff_title += ";;N/";
		diff_title += float2string(binsize_summary_histos);
		diff_title += " GeV;";

		TH1* difference;
		difference = new TH1F(diff_name.c_str(), diff_title.c_str(), binnumber,
				absolutelowerbound, absoluteupperbound);
		int nbins;
		nbins = plot_components[IDhistoA]->GetNbinsX();

		for (int jbin = 0; jbin < nbins; jbin++) {
			float sum = (plot_components[IDhistoA]->GetBinContent(jbin))
					+ (plot_components[IDhistoB]->GetBinContent(jbin));
			float diff = (plot_components[IDhistoB]->GetBinContent(jbin))
					- (plot_components[IDhistoA]->GetBinContent(jbin));

			if (sum != 0) {

				float thisbin;
				thisbin = fabs(diff / sum);
				//			thisbin = plot_components[IDhistoA]->GetBinContent(jbin);

				difference->SetBinContent(jbin, thisbin);
			}
		}

		//plot_components[IDhistoB],
		//	1, -1);

		//puts("hello");

		difference->GetYaxis()->SetRange(-1., 1.);
		vector<TH1*> vecthisto;
		vecthisto.push_back(difference);
		string deviationname;
		deviationname = plot_componentsname + "_deviation";
		TCanvas *c44;
		c44 = histo2canvas(vecthisto, plot_title, deviationname, x_label, logy,
				names);
		PrintCanvasMultiFormat(deviationname.c_str(), c44, extension,
				difference);
	}

	for (int rhisto = 0; rhisto < plot_components_unit.size(); rhisto++) {
		float area = plot_components_unit[rhisto]->Integral();
		cout << "the integral of the " << rhisto << " histogram is " << area
				<< endl;
		plot_components_unit[rhisto]->Scale(1 / area);
		cout << "the integral of the normalized " << rhisto << " histogram is "
				<< plot_components_unit[rhisto]->Integral() << endl;
	}

	gStyle->SetCanvasColor(0);
	gStyle->SetTitleColor(0);
	gStyle->SetCanvasColor(0);

	TCanvas *c3 = new TCanvas(sum_name.c_str(), sum_name.c_str(), 1);
	c3->Range(27.5, -13.2563, 152.5, 119.306);
	c3->Divide(1, 1);
	c3->cd(1);
	c3->SetBorderSize(2);
	c3->SetFrameFillColor(kWhite);
	c3->SetFillColor(kWhite);
	c3->SetBorderMode(28);
	c3->SetGrid();

	if (logy == 1) {
		gPad->SetLogy(1);
	}
	sum->SetLineColor(kRed);
	sum->Draw("");
	bg->SetLineColor(kBlue);
	bg->Draw("same");

	TLegend* legSB;
	legSB = new TLegend(0.1, 0.9 - 0.05 * 2, 0.25, 0.9);

	legSB->AddEntry(bg, "BG", "l");
	legSB->AddEntry(sum, "Sig+BG", "l");

	legSB->Draw();

	if (whichplots.sum == true) {
		PrintCanvasMultiFormat(sum_name.c_str(), c3, extension, sum);

	}

	if (whichplots.bg == true) {
		PrintCanvasMultiFormat(bg_name.c_str(), c3, extension, bg);

	}

	if (whichplots.components == true) {

		plot_componentsname += "_components";
		plot_componentsname += filenamesuffix;

		TCanvas *c22;
		c22 = histo2canvas(plot_components, plot_title, plot_componentsname,
				x_label, logy, names);
		PrintCanvasMultiFormat(plot_componentsname.c_str(), c22, extension,
				plot_components[ID_highest_histogram(plot_components)]);
	}

	if (whichplots.normalized == true) {

		plot_componentsunitname += "_components";
		plot_componentsunitname += filenamesuffix;

		TCanvas *c44;
		c44 = histo2canvas(plot_components_unit, plot_title,
				plot_componentsunitname, x_label, logy, names);
		PrintCanvasMultiFormat(
				plot_componentsunitname.c_str(),
				c44,
				extension,
				plot_components_unit[ID_highest_histogram(plot_components_unit)]);
	}
	//	PrintCanvasMultiFormat(plot_componentsunitname.c_str(), c4, extension,
	//			plot_components_unit[highest_normalized_histo_id]);

} //plothistoeps

// ************************************************
vector<vector<int> > filter_kpletlist(
		vector<vector<int> > thiscut_allONEplets_simple,
		vector<vector<int> > thiscut_classmembers,
		vector<string> single_particle_filters,
		vector<TRootLHEFParticle*> partvector,
		vector<TLorentzVector> vec_partvector) {
	// ************************************************
	//menu:
	// 1) find the filter arguments
	// 2) find the filter name
	// 3) check if the filter is a generic filter
	// 3a) in case load the needed observable
	// 4) apply the filter on the k-plets

	vector < vector<int> > last_kpletslist_kplets_filters;
	vector < vector<int> > after_filters_kpletslist;
	//cout << " found " << single_particle_filters.size() << " filters" << endl;

	last_kpletslist_kplets_filters = thiscut_allONEplets_simple;
	//cout << "class members before the filter sesction " << endl;
	//printmatrix(thiscut_classmembersPIDonly);
	for (int filter_j = 0; filter_j < single_particle_filters.size();
			filter_j++) {
		//cout << " working filter " << filter_j << endl;
		// interpret the filter string filter_name(string)
		// this amounts to call the filter passing the string
		// the filter code will valorize its parameters using the string

		// ************************************************
		// 1) find the filter arguments
		// ************************************************

		FUN thisobs_func;

		string arguments;
		arguments = get_function_name_and_arguments(
				single_particle_filters[filter_j]).second;

		// ************************************************
		// 2) find the filter name and load the function
		// ************************************************

		string filtername;
		filtername = get_function_name_and_arguments(
				single_particle_filters[filter_j]).first;

		void *dl_handletwo;
		//Dynamically open myself
		dl_handletwo = dlopen(NULL, RTLD_LAZY);

		if (!dl_handletwo) {
			cout << "Failed to open myself" << endl;
			exit(1);
		}
		char *thisfilter_funcname;
		thisfilter_funcname = const_cast<char *>(filtername.c_str());
		//Look up the function the user asked for
FIL 		thisfilter_func = (FIL) dlsym(dl_handletwo, thisfilter_funcname);
		if (!thisfilter_func) {
			cout << " Failed to find the desired filter " << filtername << endl;
			exit(1);
		}

		// ************************************************
		//		4) apply the filter on the k-plets or 1-plet
		// ************************************************

//		cout << "about to work on a list of kplet of dim="
//				<< last_kpletslist_kplets_filters.size() << endl;

		if (!last_kpletslist_kplets_filters.empty()) {

			vector < vector<int> > tempkpletslist;
			tempkpletslist.clear();
			tempkpletslist = last_kpletslist_kplets_filters;

			last_kpletslist_kplets_filters = thisfilter_func(tempkpletslist,
					thiscut_classmembers, vec_partvector, partvector,
					arguments);
		} else {
			// do nothing, as 		last_kpletslist_kplets_filters    is already empty and this will be the output
		}
		//		cout << "kplets list after the filter " << filter_j << endl;
		//		printmatrix(last_kpletslist_kplets_filters);
	} //loop on the set of filters

	after_filters_kpletslist = last_kpletslist_kplets_filters;
	// parachute for testing phase
	// thiscut_classmembers = thiscut_classmembersPIDonly;

	//
	// now succesfully tested with filter_pTmin!!!
	return after_filters_kpletslist;

}

//****************************************************************************
pair<vector<float>, vector<vector<int> > > find_the_kplet_and_compute_the_obs_on_them(
		string cut_obs, vector<int> cut_PID, int cut_k, string cut_policy,
		vector<TRootLHEFParticle*> partvector,
		vector<TLorentzVector> vec_partvector, void* dl_handle, int nmax,
		int minstatus, int maxstatus, float mcut, TFile *f, TTree * t,
		int entry, vector<int> & reusedevents, float alg,
		vector<string> single_particle_filters,
		vector<string> kplets_particle_filters) {
	//****************************************************************************
	// menu:
	// 1)  find the function for the desired observable
	// 2b) find the 1-kplets for the list made with the PID selection
	// 2a) filter according to single particle (1-plets) criteria (Energy, pT,...)
	// 2b) find the kplets for the filtered list
	// 2c) filter the list of the kplets according to k-plet criteria (deltaR, invariant mass, diff PIDS)
	// 3)  compute the observable

//	hello();

	pair < vector<float>, vector<vector<int> > > ValuesAndKplets;
	pair < vector<float>, vector<vector<int> > > ValuesAndONEplets;

	vector<float> cut_obs_kplet_value;

	//vector<vector<vector<int> > > cut_allkplets;

	vector < vector<int> > thiscut_classmembersPIDonly;
	thiscut_classmembersPIDonly.resize(0);

	vector < vector<int> > thiscut_classmembers;
	thiscut_classmembers.resize(0);

	vector < vector<vector<int> > > thiscut_allkplets;
	vector < vector<int> > thiscut_allkplets_simple;
	vector < vector<vector<int> > > thiscut_allONEplets;
	vector < vector<int> > thiscut_allONEplets_simple;

	vector < vector<int> > kpltets_after_oneplets_filter;
	vector < vector<int> > kpltets_after_kplets_filter;
	vector < vector<int> > classmembers_after_oneplets_filter;

	// *****************************************************************
	// 1) find the function for the computation of the observable
	// *****************************************************************
	char *thiscut_funcname;
	thiscut_funcname = const_cast<char *>(cut_obs.c_str());
	//Look up the function the user asked for
	//cout<< "will compute "<<thiscut_funcname<<endl;
FUN 	thiscut_func = (FUN) dlsym(dl_handle, thiscut_funcname);
	if (!thiscut_func) {
		cout << " Failed to find the desired function " << cut_obs << endl;
		exit(1);
	}

	// ************************************************************************
	// 2a) filter according to PID selection
	// ************************************************************************
//		cout << "filter according to the PID" << endl;

	thiscut_classmembersPIDonly = filter_trivial(cut_PID, vec_partvector,
			partvector, minstatus, maxstatus, mcut, f, t, entry, reusedevents,
			alg);

	// ************************************************
	// 2) find the 1-plets for this class
	// ************************************************
	//cout << "found N particles " << thiscut_classmembersPIDonly.size() << endl;

	thiscut_classmembers = thiscut_classmembersPIDonly;

	// if there is at least one object with the PID:
	// 1) apply the filters on the single particle
	// 2) find the kplet for the list after the filtering
	if (thiscut_classmembers.size() >= 1) {

		fill_allkplets_for_this_list(thiscut_allONEplets, thiscut_classmembers,
				1);

		ValuesAndONEplets.second = thiscut_allONEplets[0];
		thiscut_allONEplets_simple = thiscut_allONEplets[0];

		//	cout << "before the 1-plet filtering" << endl;
		//	printmatrix(thiscut_allONEplets_simple);
		// ************************************************************************
		// 2a) filter according to 1-plet criteria (eta, pT)
		// ************************************************************************

		//	cout << "1-filters about to be evaluated" << endl;

		kpltets_after_oneplets_filter = filter_kpletlist(
				thiscut_allONEplets_simple, thiscut_classmembers,
				single_particle_filters, partvector, vec_partvector);

		classmembers_after_oneplets_filter = oneplets2eventnumbering(
				kpltets_after_oneplets_filter, thiscut_classmembers);
		//cout << "after the 1-plet filtering" << endl;
		//	printmatrix(kpltets_after_oneplets_filter);
		//	cout << "classmembers after the 1-plet filtering" << endl;
		//	printmatrix(classmembers_after_oneplets_filter);
		// **********************************************************************
		// 2b) find the k-plets for the filtered list of single particles
		// **********************************************************************
		if (!classmembers_after_oneplets_filter.empty()) {
			// prevent from computing observables on zero particles
			// it follows an exception for the case "number" is computed
			if (classmembers_after_oneplets_filter.size() >= cut_k) {
				//		cout << " I have enough particles to make the kplets" << endl;
				//	cout << classmembers_after_oneplets_filter.size()
				//		<< "items found, required: " << cut_k << endl;

				fill_allkplets_for_this_list(thiscut_allkplets,
						classmembers_after_oneplets_filter, cut_k);
				// if cut_k=0 it returns the whole list as a single kplet

				//ValuesAndKplets.second = thiscut_allkplets[0];
				//if (cut_k > 0) {
				thiscut_allkplets_simple = thiscut_allkplets[0];
				//				hello();
				//}

				//			cout << "the kplet is made of ";
				//			printmatrix(thiscut_allkplets_simple);
				//			cout << endl;
			}

			/* this could be me moved after the size of classmembers has been checked, otherwise
			 * I am computing the combinatorics also in the case when it is not needed*/

			//cut_allkplets.push_back(thiscut_allkplets);
			//	if (thiscut_allkplets.empty() == true) {
			//		// ************************************************
			//		// **** EXCEPTION FOR THE OBSERVABLE "number" *****
			//		// ************************************************
			//		if (strcmp(cut_obs.c_str(), "number") == 0) {
			//			// in this case the result i have to add to the container is ZERO
			//
			//			cut_obs_kplet_value.push_back(0);
			//			puts("adding a zero before the single particle filters");
			//		}
			//	}
			// **********************************************************************
			// 3) filter the k-plets according to the list of filter kplets_particle_filters
			// **********************************************************************
			//		cout << "about to evaluate the k-plet filters" << endl;
			kpltets_after_kplets_filter = filter_kpletlist(
					thiscut_allkplets_simple,
					classmembers_after_oneplets_filter, kplets_particle_filters,
					partvector, vec_partvector);

//			cout <<" after the k-plet fileter the list of kplets is " <<endl;
//			printmatrix(kpltets_after_kplets_filter);

			ValuesAndKplets.second = kpltets_after_kplets_filter;

			// WHAT HAPPENS HERE IF ALL THE KPLET FAIL THE CUTS ON THE KPLETS?

			//	cout << "done with the evaluation of the k-plet filters" << endl;
		} //end of   if (!classmembers_after_oneplets_filter.empty() ) {
		else if (classmembers_after_oneplets_filter.empty()) {

			if (strcmp(cut_obs.c_str(), "number") == 0) {

				cut_obs_kplet_value.push_back(0);

				//	cout << "*** number obs exception at the stage of 1-filters"
				//		<< endl;
				//in this case the result i have to add to the container is ZERO
				//vector<vector<int> > vvitemp;
				//cut_obs_kplet_value.resize(0);
				//ValuesAndKplets.second = vvitemp;
				//puts("adding a zero after the single particle filters");
			} // end of the exception for number

		}
		// ************************************************
		// 3) compute the observable
		// ************************************************
		if (classmembers_after_oneplets_filter.empty() == false) {

			//cut_classmembers.push_back(thiscut_classmembers);
			//cout<<"particles member of the class: ";
			//printvector(cut_classmembers[icut]);
			//	puts("sono dentro il calcolo");

			if (classmembers_after_oneplets_filter.size() >= cut_k) {
				//cout << classmembers_after_oneplets_filter.size()
				//		<< "items found, required: " << cut_k << endl;
				//puts("sto per calcolare l'osservabile");
				//cut_obs_kplet_value = computeOBSonthekplet(partvector, vec_partvector, thiscut_allkplets, cut_k, thiscut_classmembers, thiscut_func);

				//	cout << "the observalble is " << endl;

				//	cout << "obesrvable " << cut_obs.c_str() << " is about to be evaluated"	<< endl;

				cut_obs_kplet_value = computeOBSonthekpletsimple(partvector,
						vec_partvector, kpltets_after_kplets_filter,
						classmembers_after_oneplets_filter, thiscut_func);

				//printvector(cut_obs_kplet_value);

			} else if (classmembers_after_oneplets_filter.size() < cut_k) {
				// ****************************************************************************************
				// The number of element found for this class are too few to compute the desired obsevable
				// Therefore I resize the result to zero. In this way the result will be an empty vector.
				// This will trigger a decision on the cut result depending on the choosen cut policy.
				// ****************************************************************************************

				//			cout << "ho trovato meno di " << cut_k << " particelle" << endl;
				cut_obs_kplet_value.resize(0);

			}

		} // end of the if  classmembers_after_oneplets_filter.empty() == false
	} // (thiscut_classmembers.size() >= 1)
	else if (thiscut_classmembers.size() < 1) {
		// no particles with the right PID have been found
		// I put an exception for the observable number

		// ************************************************
		// **** EXCEPTION FOR THE OBSERVABLE "number" *****
		// ************************************************
		if (strcmp(cut_obs.c_str(), "number") == 0) {

			cut_obs_kplet_value.push_back(0);

			//		cout << "*** number obs exception at the stage of PID" << endl;
			//in this case the result i have to add to the container is ZERO
			//vector<vector<int> > vvitemp;
			//cut_obs_kplet_value.resize(0);
			//ValuesAndKplets.second = vvitemp;
			//puts("adding a zero after the single particle filters");
		} // end of the exception for number
	}
//	cout << "the observalble is " << endl;
//	printvector(cut_obs_kplet_value);

	ValuesAndKplets.first = cut_obs_kplet_value;
	//	hello();
	//	ValuesAndKplets.second = thiscut_allkplets[0];
	//	hello();
	return ValuesAndKplets;

}
//end of the computation of the value. the size of the vector now is either zero, or

//****************************************************************************
vector<vector<float> > compute_cut_correlation_this_event(
		vector<bool> cut_passed, float w) {
	//****************************************************************************

	vector < vector<float> > cut_correlation_this_event;
	//cout << "vettore dei risultati dei cut ";
	//printvector(cut_passed);
	for (int icut = 0; icut < cut_passed.size(); icut++) {
		vector<float> jrow;

		for (int jcut = 0; jcut < cut_passed.size(); jcut++) {

			if (cut_passed[icut] == cut_passed[jcut] == 1) {
				jrow.push_back(w);
			} else {
				jrow.push_back(-w);
			}

		}
		//cout << "riga " << icut << " della correlazione ";
		//printvector(jrow);
		cut_correlation_this_event.push_back(jrow);
	}

	return cut_correlation_this_event;
}

//****************************************************************************
vector<float> compute_chained_weight_all_cuts_this_event(
		vector<bool> cut_passed, float event_weight) {
	//****************************************************************************

	vector<float> cut_bool_weigth_this_event;
	//cout << "vettore dei risultati dei cut ";
	//printvector(cut_passed);

	bool previous_are_passed;
	previous_are_passed = true;

	for (int icut = 0; icut < cut_passed.size(); icut++) {

		if ((cut_passed[icut] == 1) and (previous_are_passed)) {

			cut_bool_weigth_this_event.push_back(event_weight);

		} else if ((cut_passed[icut] == 0) or !(previous_are_passed)) {

			cut_bool_weigth_this_event.push_back(0.);
			previous_are_passed = false;

		}

	}

	return cut_bool_weigth_this_event;

}

//****************************************************************************
vector<vector<float> > compute_cut_consequentiality_this_event(
		vector<bool> cut_passed, float w) {
	//****************************************************************************

	vector < vector<float> > cut_cons_this_event;
	//cout << "vettore dei risultati dei cut ";
	//printvector(cut_passed);
	for (int icut = 0; icut < cut_passed.size(); icut++) {
		vector<float> jrow;

		for (int jcut = 0; jcut < cut_passed.size(); jcut++) {

			if ((cut_passed[icut] == 1) and (cut_passed[jcut] == 1)) {
				jrow.push_back(w);
			} else if ((cut_passed[icut] == 1) and (cut_passed[jcut] == 0)) {
				jrow.push_back(-w);
			} else {
				jrow.push_back(0);
			}

		}
		//cout << "riga " << icut << " della correlazione ";
		//printvector(jrow);
		cut_cons_this_event.push_back(jrow);
	}

	return cut_cons_this_event;
}

//****************************************************************************
//*************************************
//****************************************************************************
//*********************************************************************************************************
int main(int argc, char **argv) {
	bool debugflag;
	debugflag = false;

	bool eventoutput;
	eventoutput = false;

	ofstream debugfile;
	string debugfilename = "debuginfo";
	debugfile.open(debugfilename.c_str());
	debugfile.close();

	ofstream txteventoutput;

	/*vector<int> a;
	 a.push_back(1);
	 a.push_back(3);
	 a.push_back(5);
	 vector<int> b;
	 b.push_back(2);
	 b.push_back(4);
	 b.push_back(6);

	 vector<vector<int> > matrix;
	 matrix.push_back(a);
	 matrix.push_back(b);

	 printvector(matrix2vector(matrix));
	 */
	//	cout<< any2closestinteger(2.50)<<endl;
	if (argc != 8) {
		puts(
				"usage: ./bbar.exe observable PID1,PID2,...,PIDn k  range outputformat1,ouputformat2,...,outputformatN filename_suffix string_plot_title");
		puts(
				"./bbar.exe deltaR 93,11 2 take_all optimal .eps _JetElectron \"The title of this Plot in LaTeX\" ");
		puts(
				"./bbar.exe deltaR 93,11 2 take_all full .eps,.pdf,.root _JetElectron \"The title of this Plot in LaTeX\" ");
		puts(
				"./bbar.exe deltaR 93,11 2 take_all override,0,3.14,0.1 .svg,.txt,.root _JetElectron \"The title of this Plot in LaTeX\" ");
		puts(
				"./bbar.exe deltaR 93,11 2 take_all fixedbin,0,3.14,6.28,12.56 .svg,.txt,.root _JetElectron \"The title of this Plot in LaTeX\" ");
		exit(1);
	}
	//char blue[] = { 0x1b, '[', 'style', ';', '3', 'color', 'm', 0 };
	//style= {7 negative, 5 blink, 4 underline, 1 bold, 0 reset}
	//
	//	char blue[] = { 0x1b, '[', '1', ';', '3', '4', 'm', 0 };
	//	char ublue[] = { 0x1b, '[', '4', ';', '3', '4', 'm', 0 };
	//	char rblue[] = { 0x1b, '[', '7', ';', '3', '4', 'm', 0 };
	//	char magenta[] = { 0x1b, '[', '1', ';', '3', '5', 'm', 0 };
	//	char red[] = { 0x1b, '[', '1', ';', '3', '1', 'm', 0 };
	//	char normal[] = { 0x1b, '[', '0', ';', '3', '9', 'm', 0 };
	//	char black[] = { 0x1b, '[', '0', ';', '3', '0', 'm', 0 };
	//	char green[] = { 0x1b, '[', '0', ';', '3', '2', 'm', 0 };
	//	char yellow[] = { 0x1b, '[', '0', ';', '3', '3', 'm', 0 };
	//	char purple[] = { 0x1b, '[', '0', ';', '3', '5', 'm', 0 };
	//	char cyan[] = { 0x1b, '[', '0', ';', '3', '6', 'm', 0 };
	//	char Lgray[] = { 0x1b, '[', '0', ';', '3', '7', 'm', 0 };
	//	char Dgray[] = { 0x1b, '[', '0', ';', '3', '8', 'm', 0 };
	//	char Bred[] = { 0x1b, '[', '1', ';', '3', '1', 'm', 0 };

	/*
	 vector<vector<string> > aprocesses;
	 vector<double> asigmas;
	 analyse_proclist(aprocesses,asigmas);
	 printvector(asigmas);
	 printmatrix(aprocesses);
	 printvector(aprocesses[0]);
	 return 0;
	 */
	/*
	 vector<vector<float> > matrix;
	 vector<float> a;
	 vector<float> b, c;
	 a.push_back(-10);
	 a.push_back(-1);
	 b.push_back(100);
	 b.push_back(-0.1);
	 c.push_back(10);
	 c.push_back(-1000);
	 matrix.push_back(a);
	 matrix.push_back(b);
	 matrix.push_back(c);
	 printmatrix(matrix);
	 printmatrix(summatrix(matrix,matrix));
	 printmatrix(constant_times_matrix(10,matrix));

	 return 0;*/

	//*********************************************************************************************************
	//****************************************************************************
	//*************************************
	map<string, float> global_parameters_map;
	map<string, float> jetdefinition_parameters_map;
	map < string, string > plot_parameters_map;

	map<string, float>::iterator it;
	map<string, float>::iterator itjet;
	map<string, string>::iterator itplot;

	string signal_policy;
	signal_policy += "signal";
	void *dl_handle;
	//Dynamically open myself
	dl_handle = dlopen(NULL, RTLD_LAZY);

	if (!dl_handle) {
		cout << "Failed to open myself" << endl;
		exit(1);
	}

	char *funcname;
	string obs_name;

	// ************** CLI INTERPREATION ****************

	obs_name += argv[1];
	funcname = argv[1];

	//Look up the function the user asked for
	FUN user_func = (FUN) dlsym(dl_handle, funcname);

	//user_func = dlsym(dl_handle, funcname);
	if (!user_func) {
		cout << " Failed to find the desired function " << argv[1] << endl;
		exit(1);
	}

	cout << argc << " arguments understood" << endl;

	vector<int> classpid;

	classpid = get_CSV(argv[2]);

	//	while (numarg < argc - 3) {
	//		classpid.push_back(atoi(argv[numarg]));
	//		numarg++;
	//	}
	cout << magenta << "class is made of ";
	printvector(classpid);
	cout << normal << endl;

	//	printvector(get_CSV(argv[2]));

	int k = atoi(argv[3]);
	string req_range = argv[4];

	if (!(check_range_request(req_range))) {
		exit(1);
	}

	string extension = argv[5];
	string filenamesuffix = argv[6];
	string plot_title = argv[7];

	//	int logy, distro;
	//	if (strcmp(argv[6], "log") == 0) {
	//		logy = 1;
	//	} else if (strcmp(argv[6], "lin") == 0) {
	//		logy = 0;
	//	} else {
	//		puts("select log or lin for the Y axis");
	//		exit(1);
	//	}
	//
	//	if (strcmp(argv[7], "distro") == 0) {
	//		distro = 1;
	//	} else if (strcmp(argv[7], "events") == 0) {
	//		distro = 0;
	//	} else {
	//		puts("select distro or events for the normalization of the bins");
	//		exit(1);
	//	}

	cout << blue << argv[1] << " will be computed on the " << k
			<< "-plet(s) and plotted in the " << req_range << " range."
			<< normal << endl;
	cout << blue << "end of plot request" << normal << endl;

	// *************END CLI INTERPRETATION ****************... one day I will be a function

	/*********** CREATE OUTPUT LHE FILE */
	/*##$@#$@#$@#$ MOVE TO THE FINAL PART!!!!!!!!!!!! @!#!@$#@#$@#*/
	//ofstream homemade_lhe;
	//homemade_lhe.open("OUPUT.LHE");
	//TH2F *OutParticlePolarization = new TH2F("OutPolarization","Outgoing Particles Polarizations",53,-26,26,6,-2,2);
	gStyle->SetOptStat("nemruoi");
	//puts("ciao");
	// Load shared library
	gSystem->Load("./libExRootAnalysis.so");
	gSystem->Load("libPhysics");
	//	gSystem->Load("makedic_c.so");
	//	gSystem->Load("makevectdic_c.so");
	//	gROOT->ProcessLine(".L vectordictionary.c+");
	//	gROOT->ProcessLine(".L vvectordictionary.c+");
	gSystem->Load("vectordictionary_c.so");
	gSystem->Load("vvectordictionary_c.so");

	//	/* ************************************** */
	//	/*  THIS GENERATES THE COMBINATORIC FILE  */
	//	/* ************************************** */
	//	TFile * fileofpairings;
	//	fileofpairings = TFile::Open("pairs_in_N.root", "recreate");
	//
	//	TTree *tp;
	//	tp = new TTree("pairs", "Tree with vectors of pairs");
	//	cout << red << "Writing a new ROOT file to store the pairings " << normal << endl;
	//
	//	// Create the branches
	//	int intero;
	//	vector<vector<int> > pairings;
	//	tp->Branch("pairings", &pairings);
	//	tp->Branch("eventnumber", &intero);
	//
	//	// Generate the pairings
	//	for (int ntobeclulstered = 2; ntobeclulstered < 1000; ntobeclulstered++) {
	//		vector<vector<vector<int> > > allkplets;
	//		findkplets(allkplets, 2, ntobeclulstered);
	//		pairings = allkplets[0];
	//		intero = ntobeclulstered;
	//		//		printmatrix(pairings);
	//		//		cout << endl;
	//		cout << "trovo le coppie in " << ntobeclulstered << endl;
	//
	//		tp->Fill();
	//		tp->Write(0, TObject::kWriteDelete, 0);
	//		fileofpairings->Write();
	//	}
	//
	//	/* ************************************** */
	//	/*  THIS FINISHES THE COMBINATORIC FILE  */
	//	/* ************************************** */

	//analyse the proclist and find how many processes will be analysed
	//for each process i need to know the sigma
	//for each process I need to count how many files are given
	vector < vector<string> > processes;
	vector<double> sigmas;
	vector<double> targetcrosssections;
	vector<float> weights;
	vector<int> Nsample;
	vector < string > names;

	vector < vector<float> > global_weights;
	//[process#][event#]
	vector<float> global_total_weight;
	//[process#]

	vector < vector<float>
			> global_weights_of_the_selected_events_of_this_process;
	//[process#][event#]

	vector < vector<vector<float> > > globalresult;
	//[process#][event#][kplet#]

	vector < vector<int> > global_cut_passed;
	//[process#][cut#]

	vector < vector<float> > global_cut_efficiency;
	//[process#][cut#]

	vector<float> global_fraction_left;
	//[process#]

	vector < vector<vector<float> > > cutcorrelations_global;
	vector < vector<vector<float> > > cutconsequence_global;
	//[process#][cut#][cut#]

	vector < string > cut_obs;
	vector < vector<int> > cut_PID;
	vector<int> cut_k;
	vector<float> cut_lower;
	vector<float> cut_upper;
	vector<int> cut_howmany;
	vector < string > cut_howmany_relation;
	vector < string > cut_policy;
	vector < string > cut_IDpolicy;
	vector < string > filtri_single_particles;
	vector < vector<string> > filtri_single_particles_cut;
	vector < string > filtri_kplets;
	vector < vector<string> > filtri_kplets_cut;
	vector < string > arrangers;
	vector < vector<string> > arrangers_cut;
	bool jetrequested;
	jetrequested = false;
	struct plotsrequired whichplots;
	//**********************************************
	analyse_proclist(processes, sigmas, targetcrosssections, names, cut_obs,
			cut_PID, cut_k, cut_lower, cut_upper, cut_howmany,
			cut_howmany_relation, cut_policy, cut_IDpolicy,
			filtri_single_particles_cut, global_parameters_map,
			jetdefinition_parameters_map, plot_parameters_map,
			filtri_single_particles, filtri_kplets, filtri_kplets_cut);
	//**********************************************
	vector<double> Nsamples = sigmas;

	cout << blue << "**********************************************" << endl;
	cout << blue << "global parameters:" << normal << endl;
	for (it = global_parameters_map.begin(); it != global_parameters_map.end();
			it++) {
		cout << (*it).first << " => " << (*it).second << endl;
	}
	cout << blue << "**********************************************" << normal
			<< endl;

	cout << cyan << "**********************************************" << endl;
	cout << cyan << "plot parameters:" << normal << endl;
	for (itplot = plot_parameters_map.begin();
			itplot != plot_parameters_map.end(); itplot++) {
		cout << (*itplot).first << " => " << (*itplot).second << endl;
	}
	cout << cyan << "**********************************************" << normal
			<< endl;

	cout << red << "**********************************************" << endl;
	cout << red << "jet parameters:" << normal << endl;
	for (itjet = jetdefinition_parameters_map.begin();
			itjet != jetdefinition_parameters_map.end(); itjet++) {
		cout << (*itjet).first << " => " << (*itjet).second << endl;
	}
	cout << red << "**********************************************" << normal
			<< endl;

	cout << green << "**********************************************" << endl;
	cout << green << "single particle filters:" << normal << endl;
	printvectornl(filtri_single_particles);

	cout << endl << green << "k-plets filters:" << normal << endl;
	printvectornl(filtri_kplets);
	cout << green << "**********************************************" << normal
			<< endl;

	vector < string > x_label;
	for (int cc = 0; cc < cut_obs.size(); cc++) {
		string xx_label;
		/*cutoutput += "new ";
		 cutoutput += cut_policy[cc];
		 cutoutput += " cut ";
		 cutoutput += float2string(cut_lower[cc]);
		 cutoutput += "<=";
		 cutoutput += cut_obs[cc];
		 cutoutput += "<=";
		 cutoutput += float2string(cut_upper[cc]);
		 cutoutput += " for N";
		 cutoutput += cut_howmany_relation[cc];
		 cutoutput += int2string(cut_howmany[cc]);
		 cutoutput += " ";
		 cutoutput += int2string(cut_k[cc]);
		 cutoutput += "-plet(s) built out of PIDs: ";
		 cutoutput += vector2string(cut_PID[cc]);
		 cout << cutoutput << endl; */

		xx_label += cut_policy[cc];
		xx_label += " cut requires ";
		xx_label += xx_label += float2string(cut_lower[cc]);
		xx_label += "<=";
		xx_label += cut_obs[cc];
		xx_label += "<=";
		xx_label += float2string(cut_upper[cc]);
		xx_label += " \n ";
		//xx_label += " ";

		string cutoutput;
		cutoutput = " *) ";
		cutoutput += cut_policy[cc];
		cutoutput += "ly requires events with ";
		cutoutput += float2string(cut_lower[cc]);
		cutoutput += "<=";
		cutoutput += cut_obs[cc];
		cutoutput += "<=";
		cutoutput += float2string(cut_upper[cc]);
		//	x_label.push_back(cutoutput);
		cutoutput += " for ";
		cutoutput += cut_IDpolicy[cc];
		cutoutput += " N";
		cutoutput += cut_howmany_relation[cc];
		cutoutput += int2string(cut_howmany[cc]);
		cutoutput += " pippo group(s) of ";
		cutoutput += int2string(cut_k[cc]);
		cutoutput += " object(s) with PIDs: ";
		cutoutput += vector2string(cut_PID[cc]);

		x_label.push_back(cutoutput);
	}

	/*struct structure_name {
	 member_type1 member_name1;
	 member_type2 member_name2;
	 member_type3 member_name3;
	 .
	 .
	 } object_names;*/

	int nmaxev;
	int nbegin;
	int maxstatus;
	int minstatus;
	float lumi;
	float mcut;
	float alg;

	float printevery;
	printevery = 50000.;
	float printpassedevery;
	printpassedevery = 1000.;

	alg = 0;
	lumi = 100.0; //in 1/fb
	minstatus = 0;
	maxstatus = 0;
	nbegin = 0;
	nmaxev = 1;

	//	cout << magenta << "alog"<< alg<<endl;

	valorize(lumi, "lumi", global_parameters_map);
	valorize(nmaxev, "nmaxev", global_parameters_map);
	valorize(nbegin, "nbegin", global_parameters_map);
	valorize(minstatus, "minstatus", global_parameters_map);
	valorize(maxstatus, "maxstatus", global_parameters_map);

	valorize(printpassedevery, "printpassedevery", global_parameters_map);
	valorize(printevery, "printevery", global_parameters_map);

	valorize(mcut, "mcut", jetdefinition_parameters_map);
	valorize(alg, "algorithm", jetdefinition_parameters_map);

	string scale, plottype, sumplot, bgplot, normalizedplot, componentsplot,
			diffplot;

	string debugbool;

	valorizeS(scale, "scale", plot_parameters_map);
	valorizeS(plottype, "type", plot_parameters_map);
	valorizeS(bgplot, "bg", plot_parameters_map);
	valorizeS(sumplot, "sum", plot_parameters_map);
	valorizeS(normalizedplot, "normalized", plot_parameters_map);
	valorizeS(componentsplot, "components", plot_parameters_map);
	valorizeS(diffplot, "diff", plot_parameters_map);
	valorizeS(debugbool, "debugflag", plot_parameters_map);

	if (strcmp(debugbool.c_str(), "true") == 0) {
		debugflag = true;
	}

	//	struct plotsrequired whichplots;
	whichplots.bg = false;
	whichplots.sum = false;
	whichplots.normalized = false;
	whichplots.components = false;
	whichplots.diff = false;

	//bool dosumplot, dobgplot, donormalizedplot, docomponentsplot;
	if (strcmp(bgplot.c_str(), "true") == 0) {
		whichplots.bg = true;
	}

	if (strcmp(diffplot.c_str(), "true") == 0) {
		whichplots.diff = true;
	}

	if (strcmp(sumplot.c_str(), "true") == 0) {
		whichplots.sum = true;
	}
	if (strcmp(normalizedplot.c_str(), "true") == 0) {
		whichplots.normalized = true;
	}
	if (strcmp(componentsplot.c_str(), "true") == 0) {
		whichplots.components = true;
	}

	int logy, distro;
	if (strcmp(scale.c_str(), "log") == 0) {
		logy = 1;
	} else if (strcmp(scale.c_str(), "lin") == 0) {
		logy = 0;
	} else {
		puts("select log or lin for the Y axis");
		exit(1);
	}

	if (strcmp(plottype.c_str(), "distro") == 0) {
		distro = 1;
	} else if (strcmp(plottype.c_str(), "events") == 0) {
		distro = 0;
	} else {
		puts("select distro or events for the normalization of the bins");
		exit(1);
	}
	double sumofw = 0;
	//loop on the processes previously found and build the chains
	//then do the cuts and the analysis

	/* CHECK IF SOMEONE NEED THE JETS */
	for (int kcut = 0; kcut < cut_PID.size(); kcut++) {
		if (isinclass(cut_PID[kcut], 93) == true) {
			jetrequested = true;
		}

	}
	if (isinclass(classpid, 93) == true) {
		jetrequested = true;
	}

	cout << magenta << processes.size() << " processes understood" << normal
			<< endl;
	for (int jproc = 0; jproc < processes.size(); jproc++) {

		bool jetfile_exists;
		char jetfilename[100] = "";

		/* ****************** */
		TChain chain("LHEF");
		/* ****************** */

		//TChain *chain = new TChain("LHEF");
		//if(!FillChain(chain, argv[1])) cout<<"bad file list"<<endl;
		//chain.Add("in.root");
		/* ****************** */
		for (int kfile = 0; kfile < processes[jproc].size(); kfile++) {

			chain.Add(processes[jproc][kfile].c_str());
			cout << "Adding ROOT file " << processes[jproc][kfile] << endl;

		}
		/* ****************** */

		TFile *f;
		//		TFile *friid;
		TTree *t;

		if (jetrequested == true) {

			if (processes[jproc].size() == 1) {

				string jetfile;
				jetfile = processes[jproc][0];

				jetfile.erase(jetfile.end() - 4, jetfile.end());

				strcat(jetfilename, jetfile.c_str());

				if (alg == 0) {
					strcat(jetfilename, "JADE");
				} else if (alg == 1) {
					strcat(jetfilename, "KT");

				}

				char mcutchar[5] = "";
				sprintf(mcutchar, "%.0f", mcut);

				strcat(jetfilename, mcutchar);
				strcat(jetfilename, "jets.root");

				cout << "checking the existence of Clustering ROOT file "
						<< jetfilename << endl;
				FILE *fp = fopen(jetfilename, "r");
				if (fp) {
					// exists
					jetfile_exists = true;

					//				cout << "Found Clustering ROOT file " << jetfilename << endl;

					fclose(fp);
				} else {
					// doesnt exist
					//				cout << "NOT Found Clustering ROOT file " << jetfilename << endl;
					jetfile_exists = false;
				}

			} else {
				puts(
						"chained root files forces to redo the whole jet clustering on the fly");
				jetfile_exists = false;
			}

			if (jetfile_exists == false) {

				//			friid = TFile::Open(jetfilename, "READ");

				// (string jetfilename,string treename,vector<string> branchnames,)

				f = TFile::Open(jetfilename, "recreate");

				//			cout<<"/* creo il TTre e la struttura necessaria*/"<<endl;
				//			TTree *t;
				t = new TTree("tvec", "Tree with vectors");
				cout << red << "Writing a new " << jetfilename
						<< " ROOT file to store the clustering " << normal
						<< endl;
				//Create the branches
				int intero;
				intero = -2;
				vector<int> vint;
				vint.resize(0);
				vint.push_back(-3);

				vector < vector<int> > vvint;
				vvint.push_back(vint);

				t->Branch("clusters", &vvint);
				t->Branch("eventnumber", &intero);

				//			cout<< "ho aperto un file con il numero "<<intero<<endl;
				//			t->Write();
				t->Fill();
				t->Write(0, TObject::kWriteDelete, 0);
				//			f->Write();
				f->Write(0, TObject::kWriteDelete, 0);

			}

			//		else  {

			/* ************************************************************************** */
			/* at this point a file either existed from the beginning or has been created */
			/* ************************************************************************** */

			//			friid = TFile::Open(jetfilename, "READ");
			//						f->GetObject("tvec", t);
			f = TFile::Open(jetfilename, "update");
			cout << "a TFile has been open in update mode" << endl;

			t = (TTree*) f->Get("tvec");
			cout << "a TTree called `tvec` has been open" << endl;

			//		}

		}

		vector<int> reusedevents;

		/* ****************** */
		ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
		/* ****************** */

		Long64_t numberOfEntries = treeReader->GetEntries();
		Long64_t limit = nbegin + nmaxev;

		if (limit > 0) {
			cout << magenta << "Analisys will start at the " << nbegin + 1
					<< "th event, the last analysed event will be the " << limit
					<< "th event instead of " << numberOfEntries << "th."
					<< normal << endl;
			numberOfEntries = limit;
		}

		double w;

		//		w = lumi * 1000 * sigmas[jproc] / numberOfEntries;//sigma is assumed to be in pb

		if (targetcrosssections[jproc] != 0) {
			//when the target-crosssection is give

			w = lumi * 1000 * targetcrosssections[jproc] / numberOfEntries;
			cout << red << "EVENTS OF WEGITH 1 ARE ASSUMED (LHCO default)"
					<< normal << endl;
			cout << "rescaling to target cross-section: " << w << endl;
		} else {
			// when the rescaling is given

			w = lumi * 1000 / Nsamples[jproc];
			cout << "rescaling by weigth: " << w << endl;
		}
		// in the program I will read the event weight that is given in pb
		// for weighted events sigma is a meaningless name it is used as the number of samples
		// such that having N weighted events samples the total samples has the weight of the corresponding process

		weights.push_back(w);
		// this weights vector has one component per events sample
		/* ****************** */
		TClonesArray *branchParticle = treeReader->UseBranch("Particle");
		TClonesArray *branchEvent = treeReader->UseBranch("Event");
		/* ****************** */
		Int_t nmax;

		vector<float> weights_for_this_process;
		vector<float> weights_of_the_selected_events_of_this_process;

		float totalweight_this_process = 0.;

		vector < vector<float> > results_inthisprocess;
		vector<int> cut_passed_this_process;

		vector<float> efficiencies_this_process;
		vector<float> weights_for_cut_efficiencies_this_process;
		vector<float> weighted_efficiencies_this_process;
		vector<float> chained_weights_this_process;
		vector<float> chained_efficiencies_this_process;

		vector < vector<float> > correlations_this_process;
		vector < vector<float> > consequence_this_process;
		consequence_this_process = initialize_matrix(0, cut_obs.size(),
				cut_obs.size());
		correlations_this_process = initialize_matrix(0, cut_obs.size(),
				cut_obs.size());

		efficiencies_this_process = initialize_vector(0, cut_obs.size());
		weights_for_cut_efficiencies_this_process = initialize_vector(0,
				cut_obs.size());
		chained_weights_this_process = initialize_vector(0, cut_obs.size());

		float weight_left_this_process = 0;
		float weighted_events_left = 0;

		//cout << "initial correlation matrix for the process" << endl;
		//printmatrix(consequence_this_process);

		ofstream weightfile;
		weightfile.open("./weightfile");

		string eventobsfilename;
		eventobsfilename = "EventObservables_" + names[jproc];

		if (eventoutput) {
			txteventoutput.open(eventobsfilename.c_str());
			txteventoutput.close();
		}

		vector<float> eventsobervable;

		// *****************************************************
		// ******  LOOP ON EVENTS ******************************
		// *****************************************************

		time_t rawtime;
		struct tm * timeinfo;

		time_t before;
		time(&before);

//#pragma omp parallel for
		for (Int_t entry = nbegin; entry < numberOfEntries; ++entry) {
			Int_t entryplusone;
			entryplusone = entry + 1;
			if (is_integer(entryplusone / printevery)) {
				time(&rawtime);
				time_t now;
				timeinfo = localtime(&rawtime);
				printf("Event#: %i and the current local time and date is %s",
						entry, asctime(timeinfo));
				double dif;
				time(&now);
				dif = difftime(now, before) / 60.;
				time(&before);
				printf("It took %.2lf minutes to analyse %f events\n", dif,
						printevery);
			}

			if (debugflag == true) {
				debugfile.open(debugfilename.c_str(), ios::app);
				debugfile << "*** Event #" << entry << " ***" << endl;
				debugfile.close();
			}

			/* ****************** */
			treeReader->ReadEntry(entry);
			/* ****************** */
			if (branchParticle->GetEntries() > 0) {
				TRootLHEFEvent * thisLHEevent;
				thisLHEevent = (TRootLHEFEvent *) branchEvent->At(0);

				// in pb
				float eventweight;
				eventweight = thisLHEevent->Weight;
				// no weight attached so far to the observable "scale of PDF". results will be valid only for unweighted events
				float eventscale;
				eventscale = thisLHEevent->ScalePDF;
				eventsobervable.push_back(eventscale);
				//
				if (eventweight < 0) {
					cout << "found a negative weight event " << eventweight
							<< " at event " << entry << endl;
					exit(1);
				}
				weights_for_this_process.push_back(eventweight);
				totalweight_this_process += eventweight;

				//				cout<< " the weight in the LHE of this event is "<< eventweight<<endl;
				// WHY SHOULD I WRITE THE WEIGHTS ON A FILE?
				//weightfile << eventweight << endl;
				/* ****************** */
				nmax = branchParticle->GetEntries(); //variabile globale dell'evento
				/* ****************** */

				if (debugflag == true) {
					debugfile.open(debugfilename.c_str(), ios::app);
					debugfile << "*** Event # has " << nmax << " particles ***"
							<< endl;
					debugfile.close();
				}

				TRootLHEFParticle *part[nmax]; //globale per l'evento, voglio valorizzarla con una funzione
				TLorentzVector vec_part[nmax]; //globale per l'evento, voglio valorizzarla con una funzione

				vector<TRootLHEFParticle*> partvector; //globale per l'evento, voglio valorizzarla con una funzione
				partvector.resize(0);
				vector < TLorentzVector > vec_partvector; //globale per l'evento, voglio valorizzarla con una funzione
				vec_partvector.resize(0);

				// *****************************************************
				// ****** WRITE THE 4VECTORS AND THE LHE PARTICLES *****
				// *****************************************************
				for (Int_t icount = 0; icount < nmax; ++icount) {
					part[icount] = (TRootLHEFParticle *) branchParticle->At(
							icount);
					partvector.push_back(part[icount]);
					//		  cout<<part[icount]->Status<<endl;
					vec_part[icount].SetPxPyPzE(part[icount]->Px,
							part[icount]->Py, part[icount]->Pz,
							part[icount]->E);
					vec_partvector.push_back(vec_part[icount]);

					if (debugflag == true) {
						debugfile.open(debugfilename.c_str(), ios::app);
//						debugfile
//								<< "*** Event # has one particle with four-vector"
//								<< endl;

//						fva = fvXYZ[-124.589, 2.59137, 171.694, 215.934]
						string letter;
						letter = "";
						debugfile << "fv" << letter << "=fvXYZ["
								<< part[icount]->Px << " " << part[icount]->Py
								<< " " << part[icount]->Pz << " "
								<< part[icount]->E << "]" << endl;

						debugfile.close();
					}

				} //EoLoParticles

				//				 vector<int> tobeclustered;
				//				 tobeclustered.push_back(3);
				//				 tobeclustered.push_back(1);
				//				 tobeclustered.push_back(0);
				//				 tobeclustered.push_back(4);

				//				 printmatrix(jadejets);

				// ************************
				// ***** WORK ON CUTS *****
				// ************************
				vector<bool> cut_passed;

				for (int icut = 0; icut < cut_obs.size(); icut++) {

					// *************************************************************************
					// ******** COMPUTE THE VALUE OF THE OBSERVABLE RELEVANT FOR THE CUT *******
					// *************************************************************************

					pair < vector<float>, vector<vector<int> > > valuesANDkplets;

					vector<float> cut_obs_kplet_value;

					//					cout << "sto scrivendo esamiando il cut " << icut << endl;

					if (debugflag == true) {
						debugfile.open(debugfilename.c_str(), ios::app);

						debugfile << "cut " << icut << " about to be evaluated"
								<< endl;

						debugfile.close();
					}

					valuesANDkplets =
							find_the_kplet_and_compute_the_obs_on_them(
									cut_obs[icut], cut_PID[icut], cut_k[icut],
									cut_policy[icut], partvector,
									vec_partvector, dl_handle, nmax, minstatus,
									maxstatus, mcut, f, t, entry, reusedevents,
									alg, filtri_single_particles_cut[icut],
									filtri_kplets_cut[icut]);

					if (debugflag == true) {
						debugfile.open(debugfilename.c_str(), ios::app);

						debugfile << "cut " << icut << " has been evaluated"
								<< endl;

						debugfile.close();
					}

					// the vector will remain empty in all cases not enough class
					// members are found to compute the observable.

					cut_obs_kplet_value = valuesANDkplets.first;

					if (debugflag == true) {
						vector<float> tempcut_obs_kplet_value =
								valuesANDkplets.first;

						debugfile.open(debugfilename.c_str(), ios::app);
						debugfile << "the vector of the results for the cut #"
								<< icut << " has dimension "
								<< tempcut_obs_kplet_value.size() << "."
								<< endl;
						debugfile << "It containts the values ";
						PrintVector(debugfile, tempcut_obs_kplet_value, " ");
						debugfile << " " << endl;
						debugfile.close();
					}
					// *********************************************************
					// ********* COMPARE THE VALUE WITH THE CUT RANGE **********
					// *********************************************************

					vector<bool> kplet_pass_thecut;

					if (cut_obs_kplet_value.empty() == false) {
						for (int kpletvalue = 0;
								kpletvalue < cut_obs_kplet_value.size();
								kpletvalue++) {
							kplet_pass_thecut.push_back(
									isinrange(cut_obs_kplet_value[kpletvalue],
											cut_lower[icut], cut_upper[icut]));
						}

					}

					//					cout << "valori: ";
					//					printvector(cut_obs_kplet_value);
					if (debugflag == true) {
						debugfile.open(debugfilename.c_str(), ios::app);
						debugfile
								<< "The value is in the range specified (true=1, false=0): ";
						PrintVector(debugfile, kplet_pass_thecut, " ");
						debugfile << endl;
						debugfile.close();
					}
					// **************************************************************************************
					// ****** CHECK IF THE NUMBER OF SUCCESS IN THE CUT SATISFY THE REQUIRED RELATION *******
					// **************************************************************************************

					// the observable must have been computed at least one time

					if (cut_obs_kplet_value.empty() == false) {

						int cuthowmany;
						if (cut_howmany[icut] >= 0) {
							// at least cut_howmany[icut] particles must pass the cut
							cuthowmany = cut_howmany[icut];
							//cout << "only " << cuthowmany<< " particles must pass the cut" << endl;
						} else if (cut_howmany[icut] == -1) {

							cuthowmany = cut_obs_kplet_value.size();
							//							cout << "all the " << cuthowmany << " computed values must pass the cut" << endl;
						}
						//cout << "the relation was " << howmany(kplet_pass_thecut, 1) << cut_howmany_relation[icut] << cuthowmany << endl;

						vector < vector<int> > kpletsthatpassthecut;
						if (debugflag == true) {
							debugfile.open(debugfilename.c_str(), ios::app);

							debugfile << "The list of kplets has dimension "
									<< valuesANDkplets.second.size() << endl;
							debugfile.close();
						}

						// if there are kplets out of the selection, filter out those
						// that do not pass the cuts
						// in case that there were no kplets
						// the container will have zero dimension
						// resulting a zero output for howmany()

						if (valuesANDkplets.second.empty() == false) {

							kpletsthatpassthecut = stripthroughlist(
									valuesANDkplets.second, kplet_pass_thecut);
						}

						//						cout << "kplets: ";	printmatrix(kpletsthatpassthecut);
						//						cout<<"valori: ";printvector(cut_obs_kplet_value);
						if (pass_the_request(howmany(kplet_pass_thecut, 1),
								cut_howmany_relation[icut], cuthowmany,
								kpletsthatpassthecut, cut_IDpolicy[icut])
								== true) {

							//							cout << pass_the_request(howmany(kplet_pass_thecut, 1), cut_howmany_relation[icut], cuthowmany, kpletsthatpassthecut, cut_IDpolicy[icut]) << "=>PASSED?" << endl;
							cut_passed.push_back(1);
						} else if (pass_the_request(
								howmany(kplet_pass_thecut, 1),
								cut_howmany_relation[icut], cuthowmany,
								kpletsthatpassthecut, cut_IDpolicy[icut])
								== false) {
							//							cout << pass_the_request(howmany(kplet_pass_thecut, 1), cut_howmany_relation[icut], cuthowmany, kpletsthatpassthecut, cut_IDpolicy[icut]) << "=>NOT passed?" << endl;

							cut_passed.push_back(0);
						}

					} // end if (  cut_obs_kplet_value.empty() == false) {
					else if (cut_obs_kplet_value.empty() == true) {

						//	 puts("ho capito che non c'erano abbastanza particelle e ho messo a zero la dimensione del vettore dei risultati");

						// ***************************************
						// ********* APPLY CUT POLICY ************
						// ***************************************
						// IF THE PARTICLES TO BE SUBJECTED TO SOME CUT ARE NOT FOUND,
						// THEN THE LIST OF K-PLETS IS EMPTY AND THE LIST OF VALUES IS EMPTY AS WELL.
						// WHEN THIS HAPPEN THE CUT IS ASSUMED AS NOT PASSED OR NOT PASSSED ACCORDING
						// TO THE CUT POLICY VARIABLE
						if (cut_policy[icut].compare("strict") == 0) {
							cut_passed.push_back(0);
						} else if (cut_policy[icut].compare("loose") == 0) {
							cut_passed.push_back(1);
						}
					}

				} // for (icut)

				// *******************************************************************************
				// *************** ANALYSE CUT RESULTS TO PRODUCE THE EFFICIENCIES ***************
				// *******************************************************************************
				// ANALYZE cut_passed to produce
				// the correlation matrix entries for this event and then add it to correlations_this_process
				// consequence_this_event
				// efficiencies_this_process
				// chained_efficiencies_this_process
				// *******************************************************************************
				vector < vector<float> > correlations_this_event;
				vector < vector<float> > consequence_this_event;

				if (debugflag == true) {

					debugfile.open(debugfilename.c_str(), ios::app);
					debugfile << "Result of the cuts (not passed:" << false
							<< ", passed:" << true << ")";
					PrintVector(debugfile, cut_passed, " ");
					debugfile << endl;
					debugfile.close();
				}

				float normalization;
				normalization = 1.0 / numberOfEntries;

				// *************** Efficiencies for unweighted events ****************

				correlations_this_event = compute_cut_correlation_this_event(
						cut_passed, normalization);
				correlations_this_process = summatrix(correlations_this_process,
						correlations_this_event);

				consequence_this_event =
						compute_cut_consequentiality_this_event(cut_passed,
								normalization);
				//printmatrix(consequence_this_event);
				consequence_this_process = summatrix(consequence_this_process,
						consequence_this_event);

				//printmatrix(correlations_this_event);

				// sum cut_passed to the cumulative record of the cuts of this process cut_passed_this_process
				efficiencies_this_process = sumvector(efficiencies_this_process,
						constant_times_vector(normalization, cut_passed));
				//printvector(cut_passed);
				//printvector(efficiencies_this_process);

				if (debugflag == true) {

					debugfile.open(debugfilename.c_str(), ios::app);
					debugfile << "Weighted result of the cuts (not passed:"
							<< false << ", passed:" << true << ")";
					PrintVector(
							debugfile,
							compute_chained_weight_all_cuts_this_event(
									cut_passed, w), " ");

					debugfile << endl;
					debugfile.close();
				}

				// *************** Efficiencies for WEIGHTED EVENTS ****************

				// here I sum weights, later I shall divide them by the total cross-section to get efficiencies
				weights_for_cut_efficiencies_this_process = sumvector(
						weights_for_cut_efficiencies_this_process,
						constant_times_vector(eventweight, cut_passed));

				chained_weights_this_process = sumvector(
						chained_weights_this_process,
						compute_chained_weight_all_cuts_this_event(cut_passed,
								eventweight));

				// study this sample add cut_passed_this_process to the container

				// ****************************************
				// ****** NOW MAKE THE REQUESTED PLOT *****
				// ****************************************

				vector < vector<vector<int> > > allkplets;

				// questo classmembers qui non serve a nulla
				// vector<int> classmembers;
				// classmembers.resize(0);
				//cout << "*** CUT evaluation is OVER ***" << endl;
				//cout << "**** about to compute on the observable on this event"
				//		<< endl;

				/* ************************************************* */
				/*  OLD IMPLEMENTATION
				 * WHEN THE OBSERVABLES WAS COMPUTED IN ALL CASES, IRRESPECTIVE OF THE RESULT OF THE CUTS */
				/* ************************************************* */

//				vector<float> results_thisevent;
				/* ************************************************* */
				/*			results_thisevent = find_the_kplet_and_compute_the_obs_on_them(
				 obs_name, classpid, k, signal_policy, partvector,
				 vec_partvector, dl_handle, nmax, minstatus, maxstatus,
				 mcut, f, t, entry, reusedevents, alg,
				 filtri_single_particles, filtri_kplets).first;
				 *//* ************************************************* */
				//cout << "done computing the observable on this event" << endl;
				//
				// the computation will take a decision for the case the number of needed particles is not reached
				// using some input and will return an appropriate value.
				// it will also implement an execption for the observable "number"
				// the observable will be computed only on the kplets selected by pass the filter
				//
				// **** ADD IF PASSED ALL THE CUTS ****
				if (howmany(cut_passed, 0) == 0) {
					//result for the process under examination
					//it will change with jcount
					//cout << "risultato dei cut ";
					//printvector(cut_passed);
					//cout << "risulato calcolo osservabile ";
					//printvector(results_thisevent);
					if ((cut_obs.empty() == true)
							|| (cut_passed.empty() == false)) {
						// event is added always when no cut were requested (empty cut_obs)
						// event is not added when no cut has been computed (empty cut_passed)

						/* ########################################
						 * add a selection to choose plots or jets
						 ##########################################*/

						/*************************************************************************************************/
						/********* HAVING PASSED ALL THE CUTS NOW COMPUTE THE OBSERVABLE ***********/
						/*************************************************************************************************/

						vector<float> results_thisevent;

						/* ************************************************* */
						results_thisevent =
								find_the_kplet_and_compute_the_obs_on_them(
										obs_name, classpid, k, signal_policy,
										partvector, vec_partvector, dl_handle,
										nmax, minstatus, maxstatus, mcut, f, t,
										entry, reusedevents, alg,
										filtri_single_particles, filtri_kplets).first;
						/* ************************************************* */

						/*************************************************************************************************/
						/********* HAVING PASSED ALL THE CUTS NOW STORE THE RESULT OF THE REQUESTED OBSERVABLE ***********/
						/*************************************************************************************************/

						results_inthisprocess.push_back(results_thisevent);
						// add a unit of weight to the counter for the efficiency of the cut in this process
						weighted_events_left += 1.;
						weight_left_this_process += eventweight;
						//cout << "added " << weighted_events_left << endl;
						weights_of_the_selected_events_of_this_process.push_back(
								eventweight);
						if (eventoutput) {
							txteventoutput.open(eventobsfilename.c_str(),
									ios::app);
							txteventoutput << eventweight << " " << eventscale
									<< endl;
							txteventoutput.close();

						}
						int passed;
						passed = weighted_events_left;
						//cout << passed << endl;
						if (is_integer(
								weighted_events_left / printpassedevery)) {

							//if (true) {

							float tempeff;
							tempeff = (1. * passed) / (1. * entry);

							printf(
									"(Un)weigted Event#: %i. %i events passed the cuts, at this point the efficiency is %f",
									entry, passed, tempeff);
							cout << endl;
						}

					}
				} else if (howmany(cut_passed, 0) != 0) {
					//cout << "risultato dei cut ";
					//							printvector(cut_passed);
					//							cout << "risulato calcolo osservabile ";
					//							printvector(results_thisevent);
					//							results_inthisprocess.push_back(results_thisevent);
					//							cout << "NOT added" << endl;
				}
				//check all the cuts

				//				fill_allkplets_for_these_PIDS(nmax, classpid, part, allkplets, classmembers, minstatus, maxstatus);
				//				int nmembers = classmembers.size();
				//
				//				if (nmembers == 0) {
				//					// ***** AD-HOC SOLUTION FOR THE CASE OF "NUMBER" ****
				//					// HAS TO BE TRANSPOSED TO THE GENERAL METHOD
				//					if (strcmp(argv[1], "number") == 0) {
				//						vector<float> numberzero;
				//						numberzero.push_back(0);
				//						results_inthisprocess.push_back(numberzero);
				//						//puts("adding a zero");
				//					}
				//
				//				}
				//
				//				if (nmembers < k) {
				//					// TODO: do something according with the input for the case not enough members are found (only for cuts)
				//				}

				//				if (nmembers > 0) {
				//					// ora ho tutti gli indici delle particelle e devo calcolare gli osservabili su questi indici
				//					//le osservabili vanno prese da una lista quindi la funzione deve calcolare i valori e riempire un diagramma per un osservabile gnerico
				//					// p.es. minvof(allkplets[k=2][j]) con j da 1 alla fine di allkplets[2]
				//
				//					if (k <= nmembers) {
				//						vector<float> results_thisevent;
				//						results_thisevent = computeOBSonthekplet(part, vec_part, allkplets, k, classmembers, user_func);
				//
				//						//before adding the result to the result container
				//						//check that it passes the cuts
				//						//cout<<entry<<endl;
				//						//cout<<"contiene lo zero? "<<vectorcontains(cut_passed,0)<<endl;
				//						if (howmany(cut_passed, 0) == 0) {
				//							//result for the process under examination
				//							//it will change with jcount
				//							//cout << "risultato dei cut ";
				//							//printvector(cut_passed);
				//							//cout << "risulato calcolo osservabile ";
				//							//printvector(results_thisevent);
				//							if ((cut_obs.empty() == true) || (cut_passed.empty() == false)) {
				//								results_inthisprocess.push_back(results_thisevent);
				//								// add a unit of weight to the counter for the efficiency of the cut in this process
				//								fraction_left_this_process += normalization;
				//								//cout << "added" << endl;
				//							}
				//						} else if (howmany(cut_passed, 0) != 0) {
				//							//cout << "risultato dei cut ";
				//							//							printvector(cut_passed);
				//							//							cout << "risulato calcolo osservabile ";
				//							//							printvector(results_thisevent);
				//							//							results_inthisprocess.push_back(results_thisevent);
				//							//							cout << "NOT added" << endl;
				//
				//						}//check all the cuts

				//					}
				// }//if nmembers>0
			} //if getentries

		} //EndOf LOOP OVER EVENTS

		weighted_efficiencies_this_process = constant_times_vector(
				1. / totalweight_this_process,
				weights_for_cut_efficiencies_this_process);

		chained_efficiencies_this_process = constant_times_vector(
				1. / totalweight_this_process, chained_weights_this_process);

//		float cross_section_this_process = totalweight_this_process / Nsamples[jproc];

		float cross_section_this_process = totalweight_this_process
				* weights[jproc] / (1000 * lumi);

		/* WRITE TO FILE */
		//		f->Write();
		if (jetrequested == true) {
			cout << red << "writing the file" << normal << endl;
			t->Write(0, TObject::kWriteDelete, 0);
			//			f->Write(0, TObject::kWriteDelete, 0);
		}

		cout << blue << "### " << reusedevents.size()
				<< " clusterings have been reused ###" << normal << endl;
		cout << cyan << "the sample contains " << cross_section_this_process
				<< " pb of cross-section" << normal << endl;
		cout << "this sample contains enough events for "
				<< numberOfEntries / (cross_section_this_process * lumi * 1000)
				<< " pseduo-experiments" << endl;

		//printmatrix(correlations_this_process);
		cutcorrelations_global.push_back(correlations_this_process);
		cutconsequence_global.push_back(consequence_this_process);
		globalresult.push_back(results_inthisprocess);
		global_cut_efficiency.push_back(efficiencies_this_process);
		global_fraction_left.push_back(
				(weight_left_this_process / totalweight_this_process));
		global_weights.push_back(weights_for_this_process);
		global_total_weight.push_back(totalweight_this_process);
		global_weights_of_the_selected_events_of_this_process.push_back(
				weights_of_the_selected_events_of_this_process);
		// *********
		cout << "*** *** ***" << endl;

		//		float min_nevents = lumi * 1000 * sigmas[jproc];

		//		if (numberOfEntries < min_nevents) {
		//			cout << rblue
		//					<< "WARNING the event sample contains less events that what it should to match the desidered luminosity"
		//					<< normal << endl;
		//		}

		float event_left_this_lumi;
		event_left_this_lumi = lumi * 1000 * weight_left_this_process
				/ sigmas[jproc];

		cout << red << names[jproc] << normal << ": " << event_left_this_lumi
				<< " events left after the cuts for " << lumi << "/fb" << endl;
		cout << endl;
		float actual_number_of_events_left_inthesample = weighted_events_left;
		float relative_error_total_efficiency = 1
				/ sqrt(actual_number_of_events_left_inthesample);
		//one_digit(1.22);
		//PrintVector(std::cout, global_cut_efficiency[jproc], "  ");

		cout << "corresponding to " << actual_number_of_events_left_inthesample
				<< " +- "
				<< relative_error_total_efficiency
						* actual_number_of_events_left_inthesample
				<< " weighted events left in this sample of " << numberOfEntries
				<< endl;
		//cout << "single cut efficiencies: ";
		//printvector(global_cut_efficiency[jproc]);

		vector < vector<string> > efficiencies_table_column;
		efficiencies_table_column.push_back(uniform_length_string(cut_obs));

		efficiencies_table_column.push_back(
				vector2vectorstringprecision(
						weighted_efficiencies_this_process));
		efficiencies_table_column.push_back(
				vector2vectorstringprecision(
						chained_efficiencies_this_process));
		efficiencies_table_column.push_back(
				vector2vectorstringprecision(
						progressive_efficiencies_to_steps(
								chained_efficiencies_this_process)));
		//efficiencies_table_file.open(debugfilename.c_str(), ios::app);

		cout.unsetf(ios::floatfield);
		cout << setprecision(3);
		PrintColumns(std::cout, efficiencies_table_column);
		cout.unsetf(ios::floatfield);

//		cout << "                                              \t";
//		PrintVectorFixedChars(std::cout, cut_obs, "\t", 14);
//		cout << endl;
//
//		cout << "single cut efficiencies (using weights): \t";
//		//printvector(weighted_efficiencies_this_process);
//		PrintVectorFixedDigits(std::cout, weighted_efficiencies_this_process,
//				"\t");
//		cout << endl;
//
//		cout << "progressive total efficiency (using weights): \t";
//		//printvector(chained_efficiencies_this_process);
//		PrintVectorFixedDigits(std::cout, chained_efficiencies_this_process,
//				"\t");
//		cout << endl;
//
//		cout << "efficiency of each cut-step (using weights): \t";
//		//printvector(
//		//	progressive_efficiencies_to_steps(
//		//		chained_efficiencies_this_process));
//		PrintVectorFixedDigits(
//				std::cout,
//				progressive_efficiencies_to_steps(
//						chained_efficiencies_this_process), "\t");
//		cout << endl;

		cout.unsetf(ios::floatfield);
		cout << setprecision(3) << red << "global efficiency: "
				<< global_fraction_left[jproc] << normal << endl;
		cout.unsetf(ios::floatfield);
		//		cout << "probability of equal outcome of two cuts: " << endl;
		//		printmatrix(cutcorrelations_global[jproc]);
		//		cout << "j-th cut is immplied by the i-th: " << endl;
		//		printmatrix(cutconsequence_global[jproc]);
		cout << "*** *** ***" << endl;

	} //End Loop on the processes

	/* Done, close myself */
	dlclose(dl_handle);

	cout << endl;

	float background_events = 0;
	float signal_events = 0;

	for (int wproc = 0; wproc < global_cut_efficiency.size(); wproc++) {
		size_t found;
		found = names[wproc].find("signal");
		if (found != string::npos) {
			signal_events += lumi * 1000 * global_total_weight[wproc]
					* global_fraction_left[wproc];
		} else {
			background_events += lumi * 1000 * global_total_weight[wproc]
					* global_fraction_left[wproc];
		}
	}

	cout << " counting events the significance is "
			<< signal_events / sqrt(background_events) << endl;
	cout << " S/B is " << signal_events / background_events << endl;

	vector < vector<float> > sqrt_eff;
//[proc#][cut#]

	for (int sproc = 0; sproc < global_cut_efficiency.size(); sproc++) {
		for (int bproc = 0; bproc < global_cut_efficiency.size(); bproc++) {
			for (int ycut = 0; ycut < global_cut_efficiency[sproc].size();
					ycut++) {
				if ((sproc != bproc)
						&& (global_cut_efficiency[sproc][ycut]
								/ sqrt(global_cut_efficiency[bproc][ycut]) < 1)) {
					size_t found;
					found = names[sproc].find("signal");
					if (found != string::npos) {

						//cout << rblue << "cut " << cut_lower[ycut] << "<=" << cut_obs[ycut] << "<=" << cut_upper[ycut] << " alone lessen the significance of process " << names[sproc] << " w.r.t. process " << names[bproc] << " by a factor " << global_cut_efficiency[sproc][ycut] / sqrt(global_cut_efficiency[bproc][ycut]) << normal << endl;
					}
				}
			}
		}
	}

//plot of each sample seprately and as a whole both
//as single histogram and breaking the different pieces
//cout << globalresult.size() << endl;

	multiintelliplothisto2eps(globalresult, logy, argv, names, weights,
			global_weights,
			global_weights_of_the_selected_events_of_this_process, req_range,
			extension, x_label, distro, plot_title, filenamesuffix, whichplots);

	for (int iproc = 0; iproc < processes.size(); iproc++) {

	}
//plot of the separate samples

	return 0;

}
//main

